我正在尝试从netCDF文件中提取特定的数据集,然后将该数据转换为GeoTIFF格式。目前,我已经成功使用netCDF4提取了所需的数据,所有文件中的数据都存储为1d数组(纬度、经度、所需数据),并将它们分配给2d数组。我正在处理的netcdf文件被子集化到了一个特定的区域。然而,从这里开始我就不知道该怎么做了。通过阅读以下链接,我对geotiff转换的工作原理有一些了解:
http://adventuresindevelopment.blogspot.co.uk/2008/12/create-geotiff-with-python-and-gdal.html
这是我目前拥有的:
import netCDF4
import numpy as np
from osgeo import gdal
from osgeo import osr
#Reading in data from files and extracting said data
ncfile = netCDF4.Dataset("data.nc", 'r')
dataw = ncfile.variables["dataw"][:]
lat = ncfile.variables["Latitude"][:]
long = ncfile.variables["Longitude"][:]
n = len(dataw)
x = np.zeros((n,3), float)
x[:,0] = long[:]
x[:,1] = lat[:]
x[:,2] = dataw[:]
nx = len(long)
ny = len(lat)
xmin, ymin, xmax, ymax = [long.min(), lat.min(), long.max(), lat.max()]
xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymax, 0, -yres)
#Creates 1 raster band file
dst_ds = gdal.GetDriverByName('GTiff').Create('myGeoTIFF.tif', ny, nx, 1, gdal.GDT_Float32)
dst_ds.SetGeoTransform(geotransform) # specify coords
srs = osr.SpatialReference() # establish encoding
srs.ImportFromEPSG(3857) # WGS84 lat/long
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
dst_ds.GetRasterBand(1).WriteArray(x) # write r-band to the raster
dst_ds.FlushCache() # write to disk
dst_ds = None # save, close
上述创建geotiff的过程我主要参考了这里:
如何使用Python编写/创建GeoTIFF RGB图像文件?
我尝试了基于以下理解来写入栅格的数组:一个有3列的2D数组,其中有两个是坐标,另一个是数据。我在Snap中检查结果,得到了一张黑色的页面和沿着左侧的白线。
所以我的问题如下:
我该如何从我的netcdf文件中提取必要的地理变换数据,适当调整地理变换参数,随后使用提取的经纬度+数据w数组来写入geotiff文件?