矢栅处理(6)：点矢量提取栅格值extract

发表于 分类于 阅读次数： Waline： 本文字数： 909 阅读时长 ≈ 1 分钟

基于Python完成，(批量)用点矢量来提取栅格数据中的值。数据少也可使用Arcgis中的Extract Values to Points及multi。

思路及实现

以下部分主要转载自公众号【GIS与Climate】，该号主是我一个所的师兄，他的公众号会频繁推送他研究中GIS相关的一些处理和分析，包括深度学习，推荐关注~

思路

单个栅格数据读入内存为一个二维numpy array，提取点矢量对应的栅格值实际上是读取某经纬度对应的栅格点的值（如果是投影系则是平面坐标）。

数组array只能通过行列号索引，所以实现方法：

  • 经纬度/平面坐标 -> 栅格矩阵行列号->对应栅格值
  • 即将地理经纬度坐标变换到矩阵对应的行列号，本质上这是两个坐标空间的变换，也就是仿射变换。

案例

需预先确定数据的坐标系是否统一，需在统一参考系下才能正确运算

rasterio包中的spatial index方法提供了该功能的实现，自定义函数并通过apply实现，单个图层

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import rasterop
import geopandas as gpd
def ExtractPointValue(geometry,rs):
    '''
    geometry: geodataframe的geometry列
    rs: 栅格数据
    '''
    # geometry.xy获取点几何的经纬度
    x = geometry.xy[0][0]
    y = geometry.xy[1][0]
    row, col = rs.index(x,y)
    value = rs.read(1)[row,col]
    return value
# 数据读取及运算
rs = rasterio.open('../input/rs.tif')
gpd_points = gpd.read_file('../input/shp/test_points.shp')
gpd_points['value1'] = gpd_points['geometry'].apply(ExtractPointValue,rs=rs)