前言

前面两章部分完成了图像的基本组织和大气校正，可是对于遥感影像地理相关数据，精确的地理坐标也是非常关键的。L1级数据，可惜，并没有帮我们处理地理坐标系统。无论是打开ENVI还是用gdal工具查看图像元数据，其坐标系统都是错误的或者说是未读入的。

查阅了相关的资料和博客，尝试了一些方法。首先，使用了gdal工具直接赋值坐标系统和仿射变化。

def geoImformationProcess(inputpath,outputpath):dataset = gdal.Open(inputpath)basename = os.path.basename(inputpath)xmlpath = inputpath.replace('tiff', 'xml')DOMTree = xml.dom.minidom.parse(xmlpath)dom = DOMTree.documentElement# CenterLatitude = dom.getElementsByTagName('CenterLatitude')[0].firstChild.data # CenterLongitude = dom.getElementsByTagName('CenterLongitude')[0].firstChild.data TopLeftLatitude = eval(dom.getElementsByTagName('TopLeftLatitude')[0].firstChild.data)TopLeftLongitude = eval(dom.getElementsByTagName('TopLeftLongitude')[0].firstChild.data )BottomRightLatitude = eval(dom.getElementsByTagName('BottomRightLatitude')[0].firstChild.data )BottomRightLongitude = eval(dom.getElementsByTagName('BottomRightLongitude')[0].firstChild.data )LatDiff = abs(TopLeftLatitude - BottomRightLatitude)LonDiff = abs(TopLeftLongitude -BottomRightLongitude) ResX =   LonDiff/dataset.RasterXSizeResY =  LatDiff/dataset.RasterYSize# dataset.SetGeoTransform([TopLeftLongitude,ResX,0,TopLeftLatitude,0,-ResY])srf = osr.SpatialReference()srf.ImportFromEPSG(4326)envi_driver= gdal.GetDriverByName('Gtiff')out_ds = envi_driver.CreateCopy(outputpath, dataset,1)out_ds.SetGeoTransform([TopLeftLongitude,ResX,0,TopLeftLatitude,0,-ResY])out_ds.SetSpatialRef(srf)out_ds.FlushCache()del out_ds,dataset

进行了以下操作：

获取图像的信息，边缘经纬度等
计算单个像元的分辨率
定义坐标系统，定义仿射变换
输出具有坐标系统的新图像

存在一些问题，

图像是歪的，仿射信息还需要旋转角度
根据经纬度差计算的单个像元的分辨率是有误差的
原始的坐标系统未知
等等

经过一段时间的查阅和学习后，找到了正确的方式：

文件里的rpb(多项式模型几何校正参数)不仅可以用来进行地形校正，也是可以用来几何校正。

RPC文件

rpc文件的基本描述说明可以参考以下网址，不再细节说明。

高分六号卫星影像预处理踩坑记录（一） - 知乎 (zhihu.com)https://zhuanlan.zhihu.com/p/360829934使用Python GDAL库对高分三号全极化SAR影像进行RPC几何校正（PolSARpro格式） | 航行学园 (voycn.com)http://www.voycn.com/article/shiyongpythonRFC 22:RPC地理参考 — GDAL 文档 (osgeo.cn)https://www.osgeo.cn/gdal/development/rfc/rfc22_rpc.htmlRPCs in GeoTIFF (maptools.org)http://geotiff.maptools.org/rpc_prop.html(12条消息) （Python）使用Gdal进行批量的影像RPC正射校正_卫少东的博客-CSDN博客_gdal正射校正https://blog.csdn.net/weishaodong/article/details/121189467通过阅读上述资料可以知道以下几点：

GF卫星L1级数据几何信息都是存储在RPC(B)文件中的，一些软件可以在打开tif图像时读入辅助文件，如果没有读入那么影像就没有坐标系统和任何几何信息。
gdal是可以利用rpc文件进行正射校正的
gdal并不能自动读取rpc文件信息，并且域名也有细微的不同，需要手动去写入rpc文件信息
TiFF格式也支持RPC文件的使用
DEM文件不是必须的，默认为水平 0

几何校正

几何校正的步骤可以分为如下：

书写一个函数来解析rpb文件
将rpb文件的内容写入gdal数据集对象中，使其可识别
使用gdal库中的wrap函数进行转换

简易代码如下：

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_dict

写入：

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,dataset

校正：

def 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

可能会出现投影文件缺失问题，参考以下地址。

(13条消息) python 开发GDAL报Cannot find proj.db错误_GIS开发者的博客-CSDN博客https://hanbo.blog.csdn.net/article/details/119681691?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2~default~CTRLIST~Rate-1.pc_relevant_default&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2~default~CTRLIST~Rate-1.pc_relevant_default&utm_relevant_index=2

# 投影文件的配置
os.environ['PROJ_LIB'] = r'D:\conda\Lib\site-packages\osgeo\data\proj'

这样就完成了几何校正。

结果对比

校正后结果

可以发现，ENVI软件可以自动读取辅助文件，但是并不会主动进行校正，需要使用其正射校正的功能才能得到正确的结果。

结语

主要步骤结束，后面章节介绍模块设计和全部代码的发布。

高分影像批处理第三回——RPC文件与几何校正相关推荐

高分影像批处理第二回——辐射定标与大气校正
辐射定标遥感影像的辐射定标就是将传感器接收的无效的DN(digital number)值转化为我们所需要的辐射亮度值或者大气表观反射率.要达到这个目的,我们需要知道传感器的增益和偏移值.高分影像的g ...
看懂卫星遥感数据RPC文件
PRC(Rational Polynomial Coefficients )文件是用来存储用于遥感数据几何校正的RPC模型的文件,目前多存储成xml文件.对于遥感数据来说地理坐标的精确度是十分重要的, ...
【计算机网络高分笔记】第三章：数据链路层
[计算机网络高分笔记]第三章:数据链路层标签(空格分隔):[计算机网络] 第三章:数据链路层第三章:数据链路层 3.1 数据链路层的功能 3.2 组帧 3.3 差错控制 3.3.1 检错编码 3. ...
批量处理|基于ENVI的国产高分影像批量融合工具
从本文开始,介绍ENVI软件常用的遥感影像处理功能,并采用C#+IDL混合编程方式,实现ENVI常用功能的批量处理,对应的ENVI采用5.3版本,其他版本不保证能正常使用,尤其是5.3以下版本. 1. ...
用bat脚本批处理多个应用或文件
对于使用电脑的小伙伴,在开机时是不是要打开很多软件或文件,特别对于上班族来说,什么QQ.微信.浏览器.乱七八糟的文件.工作特定的软件... ...反正我是要打开十多个的.那么,有没有我只要打开一个就能 ...
基于ArcGIS与高分影像进行绿地变化分析
1. 需求现在有某区域的高分影像和对应范围的土地利用现状数据,现在的耕地面积比二调大了很多,需要把没有备案的土地清查出来.简单的说就是要找出哪些绿地是新增的或者由其他用地类型转化而来的. 2.总体技 ...
解决.tiff文件转.pcd文件滤波后转回.tiff文件点的序列被打乱的问题
解决.tiff文件转.pcd文件滤波后转回.tiff文件点的序列被打乱的问题解决思路 1..tiff转点云.pcd格式,定义数据结构记录对应点的三轴坐标与对应的行列信息 .tiff存放点云信息,将点 ...
云计算Python自动化运维开发实战三、python文件类型
为什么80%的码农都做不了架构师?>>> 云计算Python自动化运维开发实战三.python文件类型导语: python常用的有3种文件类型 1. 源代码 py ...
Pandas简明教程：三、Pandas文件读写
文章目录 1.CSV文件 2.Excel的读写 3.HTML文件的读写 4.其它文件(数据)类型的简单说明 5.办公自动化问题简析本系列教程教程完整目录: Pandas支持了非常丰富的文件类型(见文 ...
Python编程基础：第三十三节文件复制Copy a File
第三十三节文件复制Copy a File 前言实践前言当我们需要将一个文件中的内容复制到另一个文件中时,就需要用到copyfile()函数,该函数一共有两个参数copyfile(src, ds ...

高分影像批处理第三回——RPC文件与几何校正

前言

RPC文件

几何校正

结果对比

结语

高分影像批处理第三回——RPC文件与几何校正相关推荐

最新文章

热门文章