【无标题】

在地理信息系统（GIS）和遥感领域，栅格数据和点数据是两种常见的数据类型。栅格数据通常用于表示连续的地理现象，如地形、气温、降雨量等，而点数据则用于表示离散的地理特征，如监测站、采样点等。将栅格数据的值提取到点数据中，是空间分析中的一个常见任务，这一过程可以应用于多种应用场景，如环境监测、资源管理和城市规划等。

本文将从基础原理出发，讲解栅格数据和点数据的基本概念、栅格值提取的常用方法，以及在实际操作中可能遇到的问题和解决方案。关注公众号GeodataAnalysis，回复20240604获取示例数据和代码。

基础原理

栅格数据：栅格数据是由网格单元（像元）组成的矩阵，每个像元都有一个或多个属性值。常见的栅格数据格式包括GeoTIFF、NetCDF、HDF等。栅格数据通常用于表示连续变化的地理现象，例如高程、温度和降水量。

点数据：点数据是由一组具有地理坐标的点组成，每个点可以包含多个属性。点数据通常用于表示离散的地理特征，例如气象站、采样点和观测点。

栅格值提取：栅格值提取是指从栅格数据中获取特定点位置的属性值。常见的提取方法包括：
- 最近邻插值：直接使用距离目标点最近的栅格像元的值。
- 双线性插值：通过目标点周围的四个像元的值进行加权平均，获得更为平滑的结果。
- 样条插值：使用更复杂的数学模型，基于周围多个像元的值进行插值，适用于需要高精度的场景。

下文以TIF数据为例（NC数据和HDF数据的方法也类似，这两种数据的读取方法可见公众号），采用最近邻插值提取方法，讲解将栅格数据的值提取到点的方法。

读取点数据

本文使用的点数据采用csv格式保存，文件的lon和lat列分别记录了点的经纬度，采用pandas读取该文件，进而使用geopandas将其转为GeoDataFrame格式，具体代码如下。

df = pd.read_csv('./stations.csv')
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(x=df['lon'], y=df['lat']), crs='epsg:4326')

读取TIF数据

使用rasterio读取栅格数据，首先查看一下栅格数据的坐标系和空间范围。

src = rio.open('./PM25.tif')
src.crs, src.transform

(CRS.from_epsg(4526), 
Affine(941.0902023971238, 0.0, 33699573.3583, 0.0, -941.0902024079875, 6797904.8344))

坐标转换

将点数据的坐标系转换为栅格数据的坐标系，点数据坐标转换计算量较小，所以一般优先转换点数据的坐标系。

if src.crs != gdf.crs:
    gdf_rep = gdf.to_crs(src.crs)
else:
    gdf_rep = gdf.copy()

剔除不在栅格范围内的点

# 获取栅格的边界范围
raster_bounds = src.bounds
# 检查每个点是否在栅格范围内
gdf_rep['within_raster'] = gdf_rep.apply(
    lambda row: (raster_bounds.left <= row.geometry.x <= raster_bounds.right) and 
                (raster_bounds.bottom <= row.geometry.y <= raster_bounds.top),
    axis=1
)
gdf_rep = gdf_rep.loc[gdf_rep['within_raster'], :]

计算点对应的栅格像元的行列数

# 获取栅格数据的左上角边界的坐标以及分辨率
raster_min_x, raster_max_y = src.transform.c, src.transform.f
raster_x_res, raster_y_res = src.transform.a, - src.transform.e
# 获取每个点对应的栅格像元的行列数
gdf_rep['x_num'] = gdf_rep.geometry.apply(lambda geom: int((geom.x - raster_min_x)/raster_x_res))
gdf_rep['y_num'] = gdf_rep.geometry.apply(lambda geom: int((raster_max_y - geom.y)/raster_y_res))

提取栅格像元的值

# 栅格数据值转为数组
band = 1
raster_array = src.read(band)
# 提取值
gdf_rep['pm25'] = gdf_rep.apply(lambda row: raster_array[row.y_num, row.x_num], axis=1)
# 提取的NoData值转为None
gdf_rep['pm25'] = gdf_rep['pm25'].apply(lambda x: None if x==src.nodata else x)

结果展示