转载：遥感图像空间索引和提取像元值

遥感影像有两组坐标系。一组是图片的行列坐标，另一组是空间坐标。遥感影像中的每一个像元都对应地球上的一个地点。使用rasterio可以方便的对两组坐标系进行操作。假设我们需要获取某一地理坐标的NDVI数值，如经纬度(-119.770163586, 36.741997032)，但是我们数据是UTM坐标系的，所以我们首先要将经纬度坐标投影为UTM，我们可以使用pyproj完成投影。

作者：雪涵
链接：https://zhuanlan.zhihu.com/p/412865295
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

with rasterio.open(localname) as src:
 # Use pyproj to convert point coordinates
 utm = pyproj.Proj(src.crs) # Pass CRS of image from rasterio
 lonlat = pyproj.Proj(init='epsg:4326')
 
 lon,lat = (-119.770163586, 36.741997032)
 east,north = pyproj.transform(lonlat, utm, lon, lat)
 
 print('Fresno NDVI\n-------')
 print(f'lon,lat=\t\t({lon:.2f},{lat:.2f})')
 print(f'easting,northing=\t({east:g},{north:g})')
 
 # What is the corresponding row and column in our image?
 row, col = src.index(east, north) # spatial --> image coordinates
  print(f'row,col=\t\t({row},{col})')
 
 # What is the NDVI?
 value = ndvi[row, col]
 print(f'ndvi=\t\t\t{value:.2f}')
 
 
 # Or if you see an interesting feature and want to know the spatial coordinates:
 row, col = 200, 450
 east, north = src.xy(row,col) # image --> spatial coordinates
 lon,lat = pyproj.transform(utm, lonlat, east, north)
 value = ndvi[row, col]
 print(f'''
Interesting Feature
-------
row,col= ({row},{col})
easting,northing= ({east:g},{north:g})
lon,lat=  ({lon:.2f},{lat:.2f})
ndvi= {value:.2f}
''')