Python根据坐标点提取对应遥感影像像元值

1 背景说明

在遥感影像和栅格数据处理中，一个常见需求是：已经有一系列外部坐标点，希望从栅格影像中提取这些坐标点对应位置的像元值，并将像元值与坐标点自身的属性信息共同参与后续计算。

这个需求是来自我在做WebGIS 相关实践中遇到的一个数据处理环节。在系统中，为了将遥感影像数据与已有的坐标点数据进行整合，需要根据坐标点的坐标从遥感影像中提取对应的像元值。

这些坐标点可以来自采样点、监测点、矢量点、站点数据或其他业务数据，并不要求保持规则网格形态。关键步骤是先将坐标点转换到影像自身的投影坐标系，再换算为影像数组中的行列号，最后根据行列号读取对应像元的 DN 值。

本文整理一个基于GDAL、NumPy和pyproj的 Python 封装类，用于完成栅格影像读取、坐标转换、行列号换算和坐标点像元值提取。代码中也保留了按固定间隔生成规则采样点的辅助方法，但这不是本文方法的必要方法仅仅是扩展辅助方法。

1.1 实现目标

本文代码主要实现以下功能：

① 读取栅格影像的基本信息。
② 从指定波段读取影像数组。
③ 将经纬度坐标转换为影像投影坐标。
④ 将投影坐标转换为影像数组行列号。
⑤ 根据行列号提取对应位置的 DN 值。
⑥ 可选：按固定间隔生成规则采样点，并转换为投影坐标或经纬度。

其具体路径如图所示：

1.2 适用场景

方法适合以下场景：

① 根据一批经纬度坐标提取栅格影像 DN 值。
② 根据投影坐标提取栅格影像 DN 值。
③ 投影坐标与经纬度坐标转换。
④ 影像像素行列号与空间坐标互转。
⑤ 将点位属性与遥感影像像元值进行关联分析。
⑥ 在 WebGIS 或灾害监测类系统中，将点位数据与栅格数据进行关联。
⑦ 可选生成规则采样点，用于影像抽样或快速预览。
⑧ 输出目标坐标系为wgs84。
⑨ 遥感影像必须为正射影像。

2 依赖环境

2.1 Python 库

代码中主要使用以下库：

fromosgeo.gdalimportDataset,Open,Bandfromnumpyimportndarray,meshgrid,arangefrompyprojimportCRS,Transformer

依赖说明如下：

2.2 数据要求

输入数据通常是带空间参考信息的栅格影像，例如：

.tif .img .dat

影像最好具备以下信息：

① 可被 GDAL 正常读取。
② 包含投影信息。
③ 包含地理仿射变换参数。
④ 至少包含一个可读取波段。
⑤ 遥感影像必须为正射影像。

注：遥感影像通过 GDAL 库读取，建议选用 GDAL 原生支持的影像格式。

3 核心代码

3.1 ImageObj 类

下面是封装后的影像处理类：

fromtypingimportAnyfromosgeo.gdalimportDataset,Open,Bandfromnumpyimportndarray,meshgrid,arangefrompyprojimportCRS,TransformerclassImageObj:image_file_src:strdataset:Dataset image_height:Any image_width:Any image_band_count:Any projection:strtransform:tupledef__init__(self,image_file_src:str):self.image_file_src=image_file_src self.dataset=Open(self.image_file_src)self.image_height=self.dataset.RasterYSize self.image_width=self.dataset.RasterXSize self.image_band_count=self.dataset.RasterCount self.projection=self.dataset.GetProjection()self.transform=self.dataset.GetGeoTransform()# region 1. 正向基础函数defimage_array(self,band_id:int):band:Band=self.dataset.GetRasterBand(band_id)image_array=band.ReadAsArray()returnband,image_arraydefarray_grid(self,array:ndarray,grid_size:int=100):grid_rows,grid_cols=array.shape x_value,y_value=meshgrid(arange(0,grid_cols,grid_size),arange(0,grid_rows,grid_size))returnx_value.flatten(),y_value.flatten()defgrid_proj(self,x_value:ndarray,y_value:ndarray):x_coords=x_value*self.transform[1]+self.transform[0]y_coords=y_value*self.transform[5]+self.transform[3]returnx_coords,y_coordsdefproj_wgs84(self,x_coords:ndarray,y_coords:ndarray):src_crs_proj=CRS.from_wkt(self.projection)target_crs_proj="epsg:4326"transformer=Transformer.from_crs(src_crs_proj,target_crs_proj)lat,lon=transformer.transform(x_coords,y_coords)deltransformerreturnlon,lat# endregion# region 2. 反向基础函数defwgs84_proj(self,lon:ndarray,lat:ndarray):src_crs_proj=CRS.from_epsg(4326)target_crs_proj=CRS.from_wkt(self.projection)transformer=Transformer.from_crs(src_crs_proj,target_crs_proj)x_coords,y_coords=transformer.transform(lat,lon)returnx_coords,y_coordsdefproj_grid(self,x_coords:ndarray,y_coords:ndarray):x_value=(x_coords-self.transform[0])//self.transform[1]y_value=(y_coords-self.transform[3])//self.transform[5]returnx_value.astype(int),y_value.astype(int)# endregiondefgrid_dn(self,array:ndarray,x_value:ndarray,y_value:ndarray):points_value_array=array[y_value,x_value]returnpoints_value_arraydefarray_proj(self,array:ndarray,grid_size:int=100):x_value,y_value=self.array_grid(array,grid_size)x_coords,y_coords=self.grid_proj(x_value,y_value)returnx_coords,y_coordsdefgrid_wgs84(self,x_value:ndarray,y_value:ndarray):x_coords,y_coords=self.grid_proj(x_value,y_value)lon,lat=self.proj_wgs84(x_coords,y_coords)returnlon,latdefarray_dn(self,array:ndarray,grid_size:int=100):x_value,y_value=self.array_grid(array,grid_size)points_value_array=self.grid_dn(array,x_value,y_value)returnpoints_value_arraydefarray_wgs84(self,array:ndarray,grid_size:int=100):x_value,y_value=self.array_grid(array,grid_size)x_coords,y_coords=self.grid_proj(x_value,y_value)lon,lat=self.proj_wgs84(x_coords,y_coords)returnlon,latdefproj_dn(self,array:ndarray,x_coords:ndarray,y_coords:ndarray):x_value,y_value=self.proj_grid(x_coords,y_coords)points_value_array=self.grid_dn(array,x_value,y_value)returnpoints_value_arraydefwgs84_dn(self,array:ndarray,lon:ndarray,lat:ndarray):x_coords,y_coords=self.wgs84_proj(lon,lat)x_value,y_value=self.proj_grid(x_coords,y_coords)points_value_array=self.grid_dn(array,x_value,y_value)returnpoints_value_arraydefwgs84_grid(self,lon:ndarray,lat:ndarray):x_coords,y_coords=self.wgs84_proj(lon,lat)x_value,y_value=self.proj_grid(x_coords,y_coords)returnx_value,y_value

4 核心概念说明

4.1 image

image表示原始栅格影像文件。

在代码中，影像通过 GDAL 打开：

self.dataset=Open(self.image_file_src)

打开后可以获取影像高度、宽度、波段数、投影信息和地理变换参数。

4.2 array

array表示从影像某个波段读取出来的二维数组。

band=self.dataset.GetRasterBand(band_id)image_array=band.ReadAsArray()

二维数组中的每一个值通常对应一个像元的 DN 值。

4.3 grid

grid在本文代码中表示影像数组中的行列号位置，也就是像元索引。

它不要求空间点必须是规则网格。只要能将外部坐标点转换为影像行列号，就可以用同一套逻辑读取像元值。

对于已有坐标点，典型流程是：

wgs84/proj -> grid -> dn

代码中也提供了一个array_grid()方法，用于在没有外部点位时，按固定间隔从影像数组中生成规则采样点：

在代码中：

x_value,y_value=meshgrid(arange(0,grid_cols,grid_size),arange(0,grid_rows,grid_size))

其中：

如果grid_size=100，表示每隔 100 个像元生成一个规则采样点。这只是辅助用法，不是提取坐标点像元值的必要条件。

4.4 proj

proj表示影像自身投影坐标系下的空间坐标。

行列号转换为投影坐标时，使用的是 GDAL 的仿射变换参数：

x_coords=x_value*self.transform[1]+self.transform[0]y_coords=y_value*self.transform[5]+self.transform[3]

其中self.transform来自：

self.transform=self.dataset.GetGeoTransform()

4.5 wgs84

wgs84表示 WGS84 经纬度坐标。

代码中使用pyproj将影像投影坐标转换为 EPSG:4326：

src_crs_proj=CRS.from_wkt(self.projection)target_crs_proj="epsg:4326"transformer=Transformer.from_crs(src_crs_proj,target_crs_proj)

4.6 dn

dn表示影像像元值，也就是在指定行列号位置读取到的影像数组值。

points_value_array=array[y_value,x_value]

5 方法说明

5.1 读取影像数组

image_array(self,band_id:int)

作用：读取指定波段，并返回波段对象和影像数组。

5.2 可选生成规则采样点

array_grid(self,array:ndarray,grid_size:int=100)

作用：根据数组大小和采样间隔生成规则采样点行列号。

返回：

x_value, y_value

如果已经有外部坐标点，例如一组经纬度点，则不需要先调用这个方法。

5.3 行列号转投影坐标

grid_proj(self,x_value:ndarray,y_value:ndarray)

作用：根据 GDAL 仿射变换参数，将像元行列号转换为投影坐标。

5.4 投影坐标转经纬度

proj_wgs84(self,x_coords:ndarray,y_coords:ndarray)

作用：将影像投影坐标转换为 WGS84 经纬度坐标。

5.5 经纬度转投影坐标

wgs84_proj(self,lon:ndarray,lat:ndarray)

作用：将 WGS84 经纬度坐标转换回影像自身投影坐标系。

5.6 投影坐标转行列号

proj_grid(self,x_coords:ndarray,y_coords:ndarray)

作用：将投影坐标转换为影像数组行列号。

5.7 根据行列号读取 DN 值

grid_dn(self,array:ndarray,x_value:ndarray,y_value:ndarray)

作用：根据行列号从影像数组中提取像元值。

6 使用示例

6.1 打开影像并读取波段

image=ImageObj("example.tif")band,array=image.image_array(1)

这里读取第 1 个波段，并得到对应的二维数组。

6.2 根据经纬度点提取 DN 值

fromnumpyimportarray lon=array([113.10,113.20,113.30])lat=array([35.10,35.20,35.30])dn_values=image.wgs84_dn(array,lon,lat)

这才是本文方法最核心的用法：给定一组经纬度坐标，直接提取这些点在影像上的像元值。这些点可以任意分布，不需要保持规则网格。

6.3 根据投影坐标提取 DN 值

如果已有坐标是影像投影坐标系下的坐标，可以直接使用：

x_coords=array([385000.0,386000.0,387000.0])y_coords=array([3885000.0,3886000.0,3887000.0])dn_values=image.proj_dn(array,x_coords,y_coords)

这样可以跳过经纬度到投影坐标的转换步骤。

6.4 可选生成规则采样点

如果没有外部坐标点，只是想从影像中按固定间隔抽取一批规则采样点，可以使用：

lon,lat=image.array_wgs84(array,grid_size=100)dn_values=image.array_dn(array,grid_size=100)

这一步会根据影像数组大小生成规则行列号采样点，再提取这些采样点对应的经纬度和 DN 值。

这个用法适合影像抽样或预览，但不是“根据已有坐标点提取像元值”的必要步骤。

7 坐标转换流程

7.1 坐标点提取 DN 值流程

对于一组已有经纬度点，提取 DN 值的主流程是：

wgs84 -> proj -> grid -> dn

含义如下：

如果输入点已经是影像投影坐标，则流程可以简化为：

proj -> grid -> dn

7.2 可选规则采样点生成流程

如果没有外部坐标点，需要从影像自身生成规则采样点，可以使用：

image -> array -> grid -> proj -> wgs84

含义如下：

这一流程用于自动生成采样点，不代表输入坐标点必须是规则网格。

7.3 反查行列号流程

如果只需要得到坐标点在影像数组中的行列号，可以使用：

wgs84 -> proj -> grid

含义如下：

8 注意事项

8.1 grid_size 只用于规则采样点

grid_size只在规则采样方法时使用。如果输入是一组已有坐标点，并使用wgs84_dn()或proj_dn()提取 DN 值，则不需要设置grid_size。

当需要自动生成规则采样点时，grid_size越小，采样点越密集；grid_size越大，采样点越稀疏。

8.2 影像投影信息

如果影像缺少投影信息，CRS.from_wkt(self.projection)可能无法正常构建坐标系。因此，使用前需要确认影像包含有效的投影定义。

9 总结

本文是一个基于GDAL、NumPy和pyproj的栅格影像坐标点像元值提取方法。它的核心思想是：将一组经纬度坐标或投影坐标转换为影像数组行列号，再根据行列号读取对应位置的 DN 值。

这些输入坐标点可以任意分布，不需要保持规则网格形态。代码中按固定间隔生成规则采样点的功能只是辅助方法，适合在没有外部点位时进行影像抽样或快速预览。

实际使用时，需要重点关注影像投影信息、行列号越界和输入点是否落在影像范围内等问题。