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依赖说明如下库作用GDAL读取栅格影像、波段、投影和仿射变换参数NumPy处理影像数组和批量坐标点pyproj完成投影坐标与 WGS84 经纬度之间的转换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_srcimage_file_src self.datasetOpen(self.image_file_src)self.image_heightself.dataset.RasterYSize self.image_widthself.dataset.RasterXSize self.image_band_countself.dataset.RasterCount self.projectionself.dataset.GetProjection()self.transformself.dataset.GetGeoTransform()# region 1. 正向基础函数defimage_array(self,band_id:int):band:Bandself.dataset.GetRasterBand(band_id)image_arrayband.ReadAsArray()returnband,image_arraydefarray_grid(self,array:ndarray,grid_size:int100):grid_rows,grid_colsarray.shape x_value,y_valuemeshgrid(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_coordsx_value*self.transform[1]self.transform[0]y_coordsy_value*self.transform[5]self.transform[3]returnx_coords,y_coordsdefproj_wgs84(self,x_coords:ndarray,y_coords:ndarray):src_crs_projCRS.from_wkt(self.projection)target_crs_projepsg:4326transformerTransformer.from_crs(src_crs_proj,target_crs_proj)lat,lontransformer.transform(x_coords,y_coords)deltransformerreturnlon,lat# endregion# region 2. 反向基础函数defwgs84_proj(self,lon:ndarray,lat:ndarray):src_crs_projCRS.from_epsg(4326)target_crs_projCRS.from_wkt(self.projection)transformerTransformer.from_crs(src_crs_proj,target_crs_proj)x_coords,y_coordstransformer.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_arrayarray[y_value,x_value]returnpoints_value_arraydefarray_proj(self,array:ndarray,grid_size:int100):x_value,y_valueself.array_grid(array,grid_size)x_coords,y_coordsself.grid_proj(x_value,y_value)returnx_coords,y_coordsdefgrid_wgs84(self,x_value:ndarray,y_value:ndarray):x_coords,y_coordsself.grid_proj(x_value,y_value)lon,latself.proj_wgs84(x_coords,y_coords)returnlon,latdefarray_dn(self,array:ndarray,grid_size:int100):x_value,y_valueself.array_grid(array,grid_size)points_value_arrayself.grid_dn(array,x_value,y_value)returnpoints_value_arraydefarray_wgs84(self,array:ndarray,grid_size:int100):x_value,y_valueself.array_grid(array,grid_size)x_coords,y_coordsself.grid_proj(x_value,y_value)lon,latself.proj_wgs84(x_coords,y_coords)returnlon,latdefproj_dn(self,array:ndarray,x_coords:ndarray,y_coords:ndarray):x_value,y_valueself.proj_grid(x_coords,y_coords)points_value_arrayself.grid_dn(array,x_value,y_value)returnpoints_value_arraydefwgs84_dn(self,array:ndarray,lon:ndarray,lat:ndarray):x_coords,y_coordsself.wgs84_proj(lon,lat)x_value,y_valueself.proj_grid(x_coords,y_coords)points_value_arrayself.grid_dn(array,x_value,y_value)returnpoints_value_arraydefwgs84_grid(self,lon:ndarray,lat:ndarray):x_coords,y_coordsself.wgs84_proj(lon,lat)x_value,y_valueself.proj_grid(x_coords,y_coords)returnx_value,y_value4 核心概念说明4.1 imageimage表示原始栅格影像文件。在代码中影像通过 GDAL 打开self.datasetOpen(self.image_file_src)打开后可以获取影像高度、宽度、波段数、投影信息和地理变换参数。4.2 arrayarray表示从影像某个波段读取出来的二维数组。bandself.dataset.GetRasterBand(band_id)image_arrayband.ReadAsArray()二维数组中的每一个值通常对应一个像元的 DN 值。4.3 gridgrid在本文代码中表示影像数组中的行列号位置也就是像元索引。它不要求空间点必须是规则网格。只要能将外部坐标点转换为影像行列号就可以用同一套逻辑读取像元值。对于已有坐标点典型流程是wgs84/proj - grid - dn代码中也提供了一个array_grid()方法用于在没有外部点位时按固定间隔从影像数组中生成规则采样点在代码中x_value,y_valuemeshgrid(arange(0,grid_cols,grid_size),arange(0,grid_rows,grid_size))其中变量含义x_value像元列号y_value像元行号grid_size采样间隔如果grid_size100表示每隔 100 个像元生成一个规则采样点。这只是辅助用法不是提取坐标点像元值的必要条件。4.4 projproj表示影像自身投影坐标系下的空间坐标。行列号转换为投影坐标时使用的是 GDAL 的仿射变换参数x_coordsx_value*self.transform[1]self.transform[0]y_coordsy_value*self.transform[5]self.transform[3]其中self.transform来自self.transformself.dataset.GetGeoTransform()4.5 wgs84wgs84表示 WGS84 经纬度坐标。代码中使用pyproj将影像投影坐标转换为 EPSG:4326src_crs_projCRS.from_wkt(self.projection)target_crs_projepsg:4326transformerTransformer.from_crs(src_crs_proj,target_crs_proj)4.6 dndn表示影像像元值也就是在指定行列号位置读取到的影像数组值。points_value_arrayarray[y_value,x_value]5 方法说明5.1 读取影像数组image_array(self,band_id:int)作用读取指定波段并返回波段对象和影像数组。5.2 可选生成规则采样点array_grid(self,array:ndarray,grid_size:int100)作用根据数组大小和采样间隔生成规则采样点行列号。返回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 打开影像并读取波段imageImageObj(example.tif)band,arrayimage.image_array(1)这里读取第 1 个波段并得到对应的二维数组。6.2 根据经纬度点提取 DN 值fromnumpyimportarray lonarray([113.10,113.20,113.30])latarray([35.10,35.20,35.30])dn_valuesimage.wgs84_dn(array,lon,lat)这才是本文方法最核心的用法给定一组经纬度坐标直接提取这些点在影像上的像元值。这些点可以任意分布不需要保持规则网格。6.3 根据投影坐标提取 DN 值如果已有坐标是影像投影坐标系下的坐标可以直接使用x_coordsarray([385000.0,386000.0,387000.0])y_coordsarray([3885000.0,3886000.0,3887000.0])dn_valuesimage.proj_dn(array,x_coords,y_coords)这样可以跳过经纬度到投影坐标的转换步骤。6.4 可选生成规则采样点如果没有外部坐标点只是想从影像中按固定间隔抽取一批规则采样点可以使用lon,latimage.array_wgs84(array,grid_size100)dn_valuesimage.array_dn(array,grid_size100)这一步会根据影像数组大小生成规则行列号采样点再提取这些采样点对应的经纬度和 DN 值。这个用法适合影像抽样或预览但不是“根据已有坐标点提取像元值”的必要步骤。7 坐标转换流程7.1 坐标点提取 DN 值流程对于一组已有经纬度点提取 DN 值的主流程是wgs84 - proj - grid - dn含义如下步骤含义wgs84 - proj经纬度转换为影像投影坐标proj - grid投影坐标转换为影像行列号grid - dn根据行列号读取像元值如果输入点已经是影像投影坐标则流程可以简化为proj - grid - dn7.2 可选规则采样点生成流程如果没有外部坐标点需要从影像自身生成规则采样点可以使用image - array - grid - proj - wgs84含义如下步骤含义image - array从影像波段读取二维数组array - grid根据采样间隔生成规则行列号采样点grid - proj将行列号转换为投影坐标proj - wgs84将投影坐标转换为经纬度这一流程用于自动生成采样点不代表输入坐标点必须是规则网格。7.3 反查行列号流程如果只需要得到坐标点在影像数组中的行列号可以使用wgs84 - proj - grid含义如下步骤含义wgs84 - proj经纬度转换为影像投影坐标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 值。这些输入坐标点可以任意分布不需要保持规则网格形态。代码中按固定间隔生成规则采样点的功能只是辅助方法适合在没有外部点位时进行影像抽样或快速预览。实际使用时需要重点关注影像投影信息、行列号越界和输入点是否落在影像范围内等问题。