本文主要围绕遥感影像的 RPC（Rational Polynomial Coefficients，有理多项式系数）信息处理及地理校正相关操作展开。通过 Python 语言，利用gdal、spectral等库实现了多个功能函数，包括从 HDR 文件中读取并提取 RPC 信息、将 HDR 文件中的 RPC 信息转换并设置到栅格数据集中以进行地理校正并输出为 GeoTIFF 格式，还涵盖了对.rpc文件的解析、向 TIFF 影像写入 RPC 域信息以及进行 RPC 校正等操作，详细介绍了各功能函数的具体实现逻辑及相关参数设置情况，为遥感影像基于 RPC 信息的处理及校正提供了一套可行的代码实现方案。
import numpy as np
from osgeo import gdal
import spectraldef read_hdr_file(hdr_file):"""使用 spectral 读取 HDR 文件，并提取其中的元数据"""hdr = spectral.io.envi.read_envi_header(hdr_file)return hdrdef extract_rpc_info(rpc_list):"""格式化 RPC 信息"""rpc_metadata = {'LINE_OFF': float(rpc_list[0]),'SAMP_OFF': float(rpc_list[1]),'LAT_OFF': float(rpc_list[2]),'LONG_OFF': float(rpc_list[3]),'HEIGHT_OFF': float(rpc_list[4]),'LINE_SCALE': float(rpc_list[5]),'SAMP_SCALE': float(rpc_list[6]),'LAT_SCALE': float(rpc_list[7]),'LONG_SCALE': float(rpc_list[8]),'HEIGHT_SCALE': float(rpc_list[9]),'LINE_NUM_COEFF': list(map(float, rpc_list[10:30])),'LINE_DEN_COEFF': list(map(float, rpc_list[30:50])),'SAMP_NUM_COEFF': list(map(float, rpc_list[50:70])),'SAMP_DEN_COEFF': list(map(float, rpc_list[70:90])),}return rpc_metadata
# keys = ['ERR_BIAS', 'ERR_RAND', 'LINE_OFF', 'SAMP_OFF', 'LAT_OFF', 'LONG_OFF',
# 'HEIGHT_OFF', 'LINE_SCALE', 'SAMP_SCALE', 'LAT_SCALE',
# 'LONG_SCALE', 'HEIGHT_SCALE', 'LINE_NUM_COEFF', 'LINE_DEN_COEFF',
# 'SAMP_NUM_COEFF', 'SAMP_DEN_COEFF']def convert_hdr_rpc_to_geotiff(raster_file, hdr_file, output_tiff):# 打开栅格数据集raster_dataset = gdal.Open(raster_file)# 使用 spectral 读取并提取 HDR 文件中的元数据metadata = read_hdr_file(hdr_file)rpc_list = metadata.get('rpc info', '')if not rpc_list:raise ValueError("HDR 文件中没有找到有效的 RPC 信息")# 提取并格式化 RPC 信息rpc_metadata = extract_rpc_info(rpc_list)# 设置 RPC 元数据for key, value in rpc_metadata.items():if isinstance(value, list):raster_dataset.SetMetadataItem(key, ' '.join(map(str, value)), "RPC")else:raster_dataset.SetMetadataItem(key, str(value), "RPC")# 创建 GeoTIFF 驱动driver = gdal.GetDriverByName("GTiff")# 使用 RPC 变换进行地理校正warp_options = gdal.WarpOptions(dstSRS='EPSG:4326', rpc=True, resampleAlg='bilinear')output_dataset = gdal.Warp(output_tiff, raster_dataset, options=warp_options)# 清理资源output_dataset = Noneraster_dataset = Noneprint(f"导出完成：{output_tiff}")if __name__ == "__main__":# 输入栅格文件和 HDR 文件路径raster_file = r'JL1GF02F_PMS2_202112141423144_rad.dat'hdr_file = r'JL1GF02F_PMS2_20211214142314_rad.hdr'output_tiff = r'JL1GF02F_PMS2_20211214142314_rad.tif'# 调用函数进行转换convert_hdr_rpc_to_geotiff(raster_file, hdr_file, output_tiff)def write_rpc_to_tiff(inputpath,ap = True,outpath = None):rpc_file = inputpath.replace('tiff', 'rpb')rpc_dict = parse_rpc_file(rpc_file)if ap:# 可修改读取dataset = gdal.Open(inputpath,gdal.GA_Update)# 向tif影像写入rpc域信息# 注意，这里虽然写入了RPC域信息，但实际影像还没有进行实际的RPC校正# 尽管有些RS/GIS能加载RPC域信息，并进行动态校正for k in rpc_dict.keys():dataset.SetMetadataItem(k, rpc_dict[k], 'RPC')dataset.FlushCache()del datasetelse:dataset = gdal.Open(inputpath,gdal.GA_Update)tiff_driver= gdal.GetDriverByName('Gtiff')out_ds = tiff_driver.CreateCopy(outpath, dataset, 1) for k in rpc_dict.keys():out_ds.SetMetadataItem(k, rpc_dict[k], 'RPC')out_ds.FlushCache()del out_ds,datasetdef rpc_correction(inputpath,corrtiff,dem_tif_file = None):## 设置rpc校正的参数# 原图像和输出影像缺失值设置为0，输出影像坐标系为WGS84(EPSG:4326), 重采样方法为双线性插值（bilinear，还有最邻近‘near’、三次卷积‘cubic’等可选)# 注意DEM的覆盖范围要比原影像的范围大，此外，DEM不能有缺失值，有缺失值会报错# 通常DEM在水域是没有值的（即缺失值的情况），因此需要将其填充设置为0，否则在RPC校正时会报错# 这里使用的DEM是填充0值后的SRTM V4.1 3秒弧度的DEM(分辨率为90m)# RPC_DEMINTERPOLATION=bilinear 表示对DEM重采样使用双线性插值算法# 如果要修改输出的坐标系，则要修改dstSRS参数值，使用该坐标系统的EPSG代码# 可以在网址https://spatialreference.org/ref/epsg/32650/ 查询得到EPSG代码write_rpc_to_tiff(inputpath,ap = True,outpath = None)if dem_tif_file is None:wo = gdal.WarpOptions(srcNodata=0, dstNodata=0, dstSRS='EPSG:4326', resampleAlg='bilinear', format='Gtiff',rpc=True, warpOptions=["INIT_DEST=NO_DATA"])wr = gdal.Warp(corrtiff, inputpath, options=wo)print("RPC_GEOcorr>>>")else:wo = gdal.WarpOptions(srcNodata=0, dstNodata=0, dstSRS='EPSG:4326', resampleAlg='bilinear', format='ENVI',rpc=True, warpOptions=["INIT_DEST=NO_DATA"],transformerOptions=["RPC_DEM=%s"%(dem_tif_file), "RPC_DEMINTERPOLATION=bilinear"]) wr = gdal.Warp(corrtiff, inputpath, options=wo) print("RPC_GEOcorr>>>")## 对于全海域的影像或者不使用DEM校正的话，可以将transformerOptions有关的RPC DEM关键字删掉## 即将上面gdal.WarpOptions注释掉，将下面的语句取消注释，无DEM时，影像范围的高程默认全为0del wr
def parse_rpc_file(rpc_file):# rpc_file:.rpc文件的绝对路径# rpc_dict：符号RPC域下的16个关键字的字典# 参考网址：http://geotiff.maptools.org/rpc_prop.html；# https://www.osgeo.cn/gdal/development/rfc/rfc22_rpc.htmlrpc_dict = {}with open(rpc_file) as f:text = f.read()# .rpc文件中的RPC关键词words = ['errBias', 'errRand', 'lineOffset', 'sampOffset', 'latOffset','longOffset', 'heightOffset', 'lineScale', 'sampScale', 'latScale','longScale', 'heightScale', 'lineNumCoef', 'lineDenCoef','sampNumCoef', 'sampDenCoef',]# GDAL库对应的RPC关键词keys = ['ERR_BIAS', 'ERR_RAND', 'LINE_OFF', 'SAMP_OFF', 'LAT_OFF', 'LONG_OFF','HEIGHT_OFF', 'LINE_SCALE', 'SAMP_SCALE', 'LAT_SCALE','LONG_SCALE', 'HEIGHT_SCALE', 'LINE_NUM_COEFF', 'LINE_DEN_COEFF','SAMP_NUM_COEFF', 'SAMP_DEN_COEFF']for old, new in zip(words, keys):text = text.replace(old, new)# 以‘;\n’作为分隔符text_list = text.split(';\n')# 删掉无用的行text_list = text_list[3:-2]#text_list[0] = text_list[0].split('\n')[1]# 去掉制表符、换行符、空格text_list = [item.strip('\t').replace('\n', '').replace(' ', '') for item in text_list]for item in text_list:# 去掉‘=’key, value = item.split('=')# 去掉多余的括号‘(’，‘)’if '(' in value:value = value.replace('(', '').replace(')', '')rpc_dict[key] = valuefor key in keys[:12]:# 为正数添加符号‘+’if not rpc_dict[key].startswith('-'):rpc_dict[key] = '+' + rpc_dict[key]# 为归一化项和误差标志添加单位if key in ['LAT_OFF', 'LONG_OFF', 'LAT_SCALE', 'LONG_SCALE']:rpc_dict[key] = rpc_dict[key] + ' degrees'if key in ['LINE_OFF', 'SAMP_OFF', 'LINE_SCALE', 'SAMP_SCALE']:rpc_dict[key] = rpc_dict[key] + ' pixels'if key in ['ERR_BIAS', 'ERR_RAND', 'HEIGHT_OFF', 'HEIGHT_SCALE']:rpc_dict[key] = rpc_dict[key] + ' meters'# 处理有理函数项for key in keys[-4:]:values = []for item in rpc_dict[key].split(','):#print(item)if not item.startswith('-'):values.append('+'+item)else:values.append(item)rpc_dict[key] = ' '.join(values)return rpc_dictrpc_correction(inputpath,corrtiff,dem_tif_file = None)