遥感影像有两组坐标系。一组是图片的行列坐标,另一组是空间坐标。遥感影像中的每一个像元都对应地球上的一个地点。使用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}
''')