2014-10-22 4 views
1

저는 GDAL을 처음 사용하고 젖 힙니다.netCDF와 shapefile을 다시 투영하고 영역을 일치시키는 법

netCDF에 저장된 래스터를 셰이프 파일과 비교하려고합니다. 셰이프 파일은 netCDF가 포함하는 영역의 하위 섹션이며 데이터 세트는 약간 다른 예상을 사용합니다. 셰이프 파일 데이터 세트를 netCDF 투영으로 변환합니다. netCDF 파일에는 래스터 배열과 lat, lon, x, y에 대한 1 차원 배열이 들어 있습니다.

지금 당장 제 코드는 gdal.RasterizeLayer를 사용하여 셰이프 파일을 tiff로 래스터 화 한 다음 gdal.ReprojectImage를 사용하여이를 새 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

나는 단지 다른 일을했지만 코드는 다른 컴퓨터에 있습니다. 1) rasterizeLayer를 호출하기 전에 래스터와 모양 파일을 동일한 투영에 배치해야하므로 래스터 화 후에 재 투영 할 필요가 없어야하며 2) rasterizeLayer는 래스터화할 피쳐를 파악할 것입니다. 3) 래스터의 데이터를 다시로드하지 않고 수십 개의 배열로 가져올 수 있습니다. 나중에이 문제를 다시 해결하려고 노력할 것입니다. – user1269942

답변

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은 래스터의 공간 위치 정보를 보유합니다).