私はGDALを初めて使用し、足を濡らしています。
netCDF に保存されているラスターをシェープファイルと比較しようとしています。シェープファイルは、netCDF によってカバーされる領域のサブセクションであり、データセットはわずかに異なる投影法を使用します。シェープファイル データセットを netCDF プロジェクションに変換します。netCDF ファイルには、緯度、経度、x、y のラスター配列と 1 次元配列が含まれています。
現在、私のコードは gdal.RasterizeLayer を使用してシェープファイルを tiff にラスタライズし、次に gdal.ReprojectImage を使用してそれを新しい tiff に再投影しています。私の問題は、netCDF データのサブセクションを選択する必要がある 2 番目の tiff の範囲を決定する方法がわからないことです。
私のコードの関連セクションは次のとおりです。
#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 を再定義する必要があると考えています。それを調べています。