如何重新投影和区域匹配 netCDF 和 shapefile答案

【问题标题】：How to reproject and area match a netCDF and shapefile如何重新投影和区域匹配 netCDF 和 shapefile
【发布时间】：2014-10-22 19:50:40
【问题描述】：

我是 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 以进行投影转换。调查一下。

【问题讨论】：

只是一个简短的评论，我前几天刚做了这个，但代码在我的另一台电脑上。 1）我认为光栅和形状文件在调用 rasterizeLayer 之前需要在同一个投影中，所以你不应该在光栅化后重新投影，2） rasterizeLayer 将找出要光栅化的特征。 3）您可以将栅格数据作为 numpy 数组获取，而无需重新加载。稍后会尝试回到这个问题。

标签： python gis raster gdal osgeo

【解决方案1】：

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

#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 保存栅格的空间位置信息）。

【讨论】：