Как перепрограммировать и область соответствует netCDF и шейп-файлу

Question

Я новичок в GDAL и просто мокрые ноги.

Я пытаюсь сравнить растры, хранящиеся в netCDF, с шейп-файлами. Формированные профили являются подразделами области, покрываемой netCDF, а на наборах данных используются несколько разные прогнозы. Я конвертирую набор данных шейп-файла в проекцию netCDF. Файл netCDF содержит растровый массив и 1d массивы для lat, lon, x, y.

Прямо сейчас, в моем коде используется gdal.RasterizeLayer, чтобы растрировать шейп-файл в 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']) 

Re-проецировании

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 и интерполировать его для соответствия.

EDIT: Теперь я думаю, что мне нужно индивидуально преобразовать экстенты и использовать это, чтобы переопределить GeoTransform для преобразования проекции. Вглядываясь в него.

Answer 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 хранит информацию о пространственном местоположении для растра).