如何使用python从x，y，z导出到geotiff答案

【问题标题】：How to export from x,y,z to geotiff using python如何使用python从x，y，z导出到geotiff
【发布时间】：2020-12-23 17:23:54
【问题描述】：

我在 python 中有代码来打开 netcdf 并导出 x,y, 浓度。我如何使用 python 从这个变量 x,y, contration 中生成 geotiff？

提前谢谢你

代码如下：

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
from numpy.random import uniform, seed
from osgeo import gdal, osr, gdal_array
import xarray as xr
import numpy as np
import netCDF4 as nc
from datetime import datetime, timedelta
import os
import sys
import scipy
from scipy import interpolate,misc

from scipy.interpolate import griddata,RegularGridInterpolator
in_filename='E:/Ioannis_MarineEO_Task1/03_Dati_CMEMS/arctic/20190130_dm-metno-MODEL-topaz4-ARC-b20190129-fv0/20190130_dm-metno-MODEL-topaz4-ARC-b20190129-fv02.0.nc'
var_name='salinity'

error_code=0
#Open netCDF file with gdal 
src_ds=gdal.Open(in_filename)


subdataset='NETCDF:"'+in_filename+'":'+var_name    


#Read data using xarray
xr_ensemble = xr.open_dataset(in_filename)
data = xr_ensemble[var_name]
print ('Properties of the variable to be exported is: '+str((data)))
data = np.ma.masked_array(data)
data=np.ma.getdata(data)
data1=np.asarray(data)

ncDataset = nc.Dataset(in_filename, 'r')
latit = ncDataset.variables["latitude"][:]

latit=np.ma.getdata(latit)
longi = ncDataset.variables["longitude"][:]

longi=np.ma.getdata(longi)

if np.ndim(data)<=3:
    data=data[0,:,:]
else:
    data=data[0,0,:,:]
var_data=[]
x=[]
y=[]
var_data=[]


for (i,j), value in np.ndenumerate(data1):
    if np.isnan(data1[i,j])==False:
        var_data.append(data1[i,j])
        x.append(longi[i,j])
        y.append(latit[i,j])


temp1=np.column_stack([x,y,var_data])

我用 xarray 打开来自 netcdf 的数据。使用 netcdf 我打开 y（纬度）和 x（经度）。然后为每个 x,y,pixel_value 创建列表。然后我想从 x,y,pixel_value 的点使用 python 创建一个 geotif？

【问题讨论】：

到目前为止有什么尝试？你面临的问题是什么？你能分享你的代码吗？
代码如下：

标签： python

【解决方案1】：

在 Python 中，使用 rasterio 和 numpy 从 netCDF4 对象制作 geotiff 非常方便。假设您的 netcdf 文件具有 2d 或 3d 数组，则可以执行以下操作。您将需要将分配给文件的 geotiff 覆盖范围、投影、分辨率等信息。

import netCDF4 as nc
import rasterio as rio
import numpy as np

fpath = 'E:/Ioannis_MarineEO_Task1/03_Dati_CMEMS/arctic/20190130_dm-metno-MODEL-topaz4-ARC-b20190129-fv0/20190130_dm-metno-MODEL-topaz4-ARC-b20190129-fv02.0.nc'

ds = nc.Dataset(fpath)
# Check file and dataset properties by printing it.
print(ds)

# Dimensions and time

for dim in ds.dimensions.values():
    print(dim)

print(ds['time'])

变量详细信息应提供文件的最小和最大纬度和经度。

# Varaible details 
for var in ds.variables.values():
    print(var)


ds_sal = ds['salinity']

# numpy array of your salinity data 
sal_array = np.array(ds_sal)

print(sal_array.shape) # This should provide information about height and 
                         #width to make geotiff.

根据您的文件属性设置来源。对于此示例，我正在考虑该文件包含具有 0.05 度分辨率的地理投影的全球数据覆盖范围。数据的高度和宽度分别为3600和7200。

# create geotiff using above array and rasterio

from rasterio.transform import Affine, from_origin

transform = from_origin(-180,90, 0.05, 0.05)

ds_tif = rio.open('netcdf_to_geotif.tiff',
                   'w',
                    driver = 'GTiff',
                    height = 3600,
                    width = 7200,
                    dtype = sal_array.dtype,
                    count = 1,
                    crs = 'EPSG:4326',
                    transform = transform)
ds_tif.write(sal_array, 1)
ds_tif.close()

根据您的要求调整范围和分辨率，这应该会提供一个 geotiff 文件。

【讨论】：