2

我是 GDAL 的新手，我的脚刚刚湿透。

我正在尝试将存储在 netCDF 中的栅格与 shapefile 进行比较。shapefile 是 netCDF 覆盖区域的子部分，数据集使用略有不同的投影。我将 shapefile 数据集转换为 netCDF 投影。netCDF 文件包含栅格数组和 lat、lon、x、y 的一维数组。

现在，我的代码使用 gdal.RasterizeLayer 将 shapefile 光栅化为 tiff，然后使用 gdal.ReprojectImage 将其重新投影到新的 tiff 中。我的问题是我无法弄清楚如何确定第二个 tiff 的范围 - 我需要选择 netCDF 数据的子部分。

这是我的代码的相关部分：

#Extract projection information
obs_driver = ogr.GetDriverByName('ESRI Shapefile')  
obs_dataset = obs_driver.Open(obsfiles[0])
layer = obs_dataset.GetLayer()
obs_proj = layer.GetSpatialRef()
mod_proj = obs_proj.SetProjParm('parameter',90)   #Hardcode param difference
xmin,xmax,ymin,ymax = layer.GetExtent()   #Extents pre-reproject

光栅化 

tiff1 = gdal.GetDriverByName('GTiff').Create('temp1.tif', 1000, 1000, 1, gdal.GDT_Byte)
tiff1.SetGeoTransform((xmin, pixelWidth, 0, ymin, 0, pixelHeight))
gdal.RasterizeLayer(tiff1,[1],layer,options = ['ATTRIBUTE=attribute'])

重新投影

dst1 = gdal.GetDriverByName('GTiff').Create('temp3.tif', 1000, 1000, 1, gdal.GDT_Byte)
dst1.SetGeoTransform((xmin, pixelWidth, 0, ymin, 0, pixelHeight))
dst1.SetProjection(mod_proj.ExportToWkt())
gdal.ReprojectImage(gdal.Open('./temp1.tif'), dst1, obs_proj.ExportToWkt(), mod_proj.ExportToWkt(), gdalconst.GRA_Bilinear)

并将栅格转换为数组以进行逐点比较

import matplotlib.pyplot as plt
obs = plt.imread('temp3.tif')

所以现在我需要获取这个数组的范围（在新投影中），这样我就可以将它与 netCDF 数组的正确子部分匹配，并对其进行插值以匹配。

编辑：现在我认为我需要单独转换范围并使用它来重新定义 GeoTransform 以进行投影转换。调查它。

4

1 回答 1

0

想通了 - 您需要在转换数据之前将范围转换为新的投影类型。

#Calculate new x,y positions of extents
tx = osr.CoordinateTransformation(obs_proj,mod_proj)
  #Projection corrected extent of area in question
corr_ext = [0,0,0,0]  #[min,max,ymin,ymax]
(corr_ext[0],corr_ext[2],ignore_z) = tx.TransformPoint(ext[0],ext[2]) #Ignore z componenet - data is 2D
(corr_ext[1],corr_ext[3],ignore_z) = tx.TransformPoint(ext[1],ext[3])

然后使用这些范围来设置包含重投影数据的对象的 GeoTransform。（GeoTransform 保存栅格的空间位置信息）。

于 2014-10-24T14:58:28.033 回答