GIS与Climate

V1

2022/08/06阅读:0主题:默认主题

GIS数据处理:用Python提取站点处的栅格值

在进行地理数据处理的时候,我们经常需要用点矢量来提取栅格数据中的值,在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')
image
image

可视化一下:

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()

可以看到已经成功提取了各个点处的值。

总结

  1. 通过点来提取栅格处的值本质就是进行矩阵元素的索引;
  2. 巧妙利用pandas的apply函数进行批量处理;
  3. 注意矢量和栅格的坐标系要一致,不一致的需要先统一。

小思考

  1. 如果点超级多怎么办?
  • 栅格运算永远比矢栅叠加运算快,先把点栅格化再获取索引就更快了
  1. 如何提取矢量多边形的值?
  • 下篇文章。
  1. 如何批量提取栅格时间序列的值?
  • 上面的函数改一下,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/

分类:

后端

标签:

数据结构与算法

作者介绍

GIS与Climate
V1

公众号:GIS与Climate,欢迎关注