HDF 파일에서 long lat 값을 가져 오는 방법은 무엇입니까?

Question

저는 gdal을 사용하여 HDF file의 하위 데이터 세트에서 위도와 경도의 값을 얻으려고합니다. 그러나 나는 다음과 같은 오류를 받고 있어요 :

IndexError: index -62399 is out of bounds for axis 1 with size 4800 

여기 내 코드입니다 :

from osgeo import ogr, osr,gdal 

hdf_file = gdal.Open("MOD13Q1.A2017321.h31v10.006.2017337222145.hdf") 
subDatasets = hdf_file.GetSubDatasets() 
val_dict = {} 
#print subDatasets[0] 
dataset = gdal.Open(subDatasets[1][0]) 
transf = dataset.GetGeoTransform() 
success,transfInv = gdal.InvGeoTransform(transf) 
ds = dataset.ReadAsArray() 
#lon,lat = -17.586972, 139.158043 
lat = -16.718853 
lon = 142.645773 
px, py = gdal.ApplyGeoTransform(transfInv, lon, lat) 

value = ds[int(px),int(py)] 
print value 

사람이 내가 잘못 뭘하는지 말해 줄래?

Answer 1

데이터 집합의 geotransform (transf)을 보면 좌표가 위도/경도 (사인 곡선)가 아니라는 것을 알 수 있습니다.

따라서 픽셀 좌표로 변환하기 위해 지오 변환을 적용 할 때 위도/경도 값을 제공해서는 안됩니다. 이 값은 데이터 세트와 동일한 프로젝션에 있어야합니다. 당신은 왼쪽 상단 모서리의 좌표를 입력하면

예를 들어, 결과로 (0,0)을 얻을 것이다 :

mapx = transf[0] + transf[1] * ds.shape[1] 
mapy = transf[3] + transf[5] * ds.shape[0] 
px, py = gdal.ApplyGeoTransform(transfInv, mapx, mapy) 

:

mapx = transf[0] 
mapy = transf[3] 
px, py = gdal.ApplyGeoTransform(transfInv, mapx, mapy) 

또는 오른쪽 하단 모서리에

(4800, 4800)이됩니다.