在进行地理数据处理的时候，我们经常需要用点矢量来提取栅格数据中的值，在Arcgis或者R里面，这个很好实现，用extract工具或者函数就行，但是我看了下貌似Python里面没有现成的。那就自己写一个吧。
本文示例用一个点矢量的shp文件来批量提取一个栅格数据上各点处的值。
公众号《GIS与Climate》，欢迎关注和交流。
读取数据
分别用rasterio和geopandas包读取栅格和矢量数据：
rs = rasterio.open('../input/rs_min_1km.tif')
shp = gdp.read_file('../input/shp/test_points_proj.shp')
可视化一下：
fig, ax = plt.subplots()
shp.plot(ax=ax,color='orangered')
show(rs,ax=ax)
shp文件通过geopandas读取进来之后是一个geodataframe格式的数据，不要多想，直接把它看做是一个dataframe格式的就行（类似于R语言中的sf包），这样子pandas dataframe上可以用的方法我们也可以直接用在geodataframe上，那么批量处理就方便了。
从栅格提取值的本质
要知道，栅格数据在本质上就是一个array，如果想要知道某个点的值，其实就是获取这个array中某一个元素的值。如果是一个简单的numpy array，那么知道行列号就可以直接通过索引来获取值了，但是对于地理数据，我们一般是根据地理坐标来提取对应位置的栅格值，而不是行列号。
知道了本质就好办了，把对应的地理坐标变换到矩阵对应的行列号就可以了，本质上这是两个坐标空间的变换，也就是仿射变换。大佬们已经写好了这个功能，直接用rasterio包中的index方法即可，具体见参考链接【1】。
自定义提取函数
然后我们根据上面的原理写一个函数：
def ExtractPointValue(geometry,rs):
'''
geometry: geodataframe的geometry列
rs: 栅格数据
'''
x = geometry.xy[0][0]
y = geometry.xy[1][0]
row, col = rs.index(x,y)
value = rs.read(1)[row,col]
return value
批量提取值
有了函数之后，我们用apply函数直接把上述函数作用于geodataframe的geometry列即可：
shp['value1'] = shp['geometry'].apply(ExtractPointValue,rs=rs)
shp.head()
可以看到已经成功提取了各个点处的值。
总结
- 通过点来提取栅格处的值本质就是进行矩阵元素的索引；
- 巧妙利用pandas的apply函数进行批量处理；
- 注意矢量和栅格的坐标系要一致，不一致的需要先统一。
小思考
- 如果点超级多怎么办？
- 栅格运算永远比矢栅叠加运算快，先把点栅格化再获取索引就更快了
- 如何提取矢量多边形的值？
- 下篇文章。
- 如何批量提取栅格时间序列的值？
- 上面的函数改一下，read函数直接读取多波段栅格就行（不加1），或者加个循环；
- 不是多波段的图像，再用apply函数，换个参数就行；
公众号《GIS与Climate》，欢迎关注和交流。
参考
【1】https://rasterio.readthedocs.io/en/latest/quickstart.html?highlight=index#spatial-indexing
【2】https://hatarilabs.com/ih-en/extract-point-value-from-a-raster-file-with-python-geopandas-and-rasterio-tutorial
【3】https://gis.stackexchange.com/questions/260304/extract-raster-values-within-shapefile-with-pygeoprocessing-or-gdal
【4】https://www.earthdatascience.org/courses/use-data-open-source-python/spatial-data-applications/lidar-remote-sensing-uncertainty/extract-data-from-raster/
本文由mdnice多平台发布
