使用Python GDAL库对高分三号全极化SAR影像进行RPC几何校正（PolSARpro格式）

前言
数据集
- 高分三号影像
- - 高分三号数据文件结构
- DEM数据
辐射定标处理
RPC几何校正
- 基本原理
- 思路
- 实现
- - 陆地区域
  - 濒海区域
  - - 填充DEM缺失值
- 注意事项
后续处理
- S2矩阵
- - PolSARpro在坐标系上的bug
  - 校正效果检验
- T3矩阵
后语
参考文献

（原创文章，码字不易，转载请注明来源，谢谢！）

前言

很久没有更新博客，整个研三来基本都在忙着各种事情（考博、写毕业论文、找工作、党建等），现在基本收尾了。在毕业答辩之前，终于零星有点空余时间写下这篇博客，希望对国产卫星的应用推广有所帮助。

尽管高分三号SAR卫星(GF-3)早在2016年8月发射升空，但长期以来，由于高分三号获取困难以及处理工具较少，其应用相对国外SAR卫星而言较为有限。回想起本科四年级时接触GF-3影像时，那时在研究极化SAR影像的分类应用，在高分三号的预处理步骤中也是碰到很多麻烦，很多时候卡在了几何校正(也可做地理编码)这个步骤。在经过近3年开源RS/GIS的技术探索，现在在借助开源的力量，终于有能力解决这个制掣之处了。相信许多想利用PolSARpro 对高分三号SAR影像进行极化处理的同行，也受到这个困扰，希望本文能给您们带来一点帮助。

使用GDAL库和RPC文件校正卫星影像，已有不少的实现代码，但基本是都是利用C\C++等语言的实现的，相对而言，需要配置GDAL的开发环境，难度较大，且要解析针对不同遥感传感器的rpc文件记录的参数也是一个较麻烦的事情。使用python的一大好处就是拥有众多便捷的包和库，并且安装和管理便捷。因此，使用python GDAL库实现卫星影像的RPC几何校正是可行的，并且还要简单许多。

尽管本文是针对高分三号卫星为例的，但是只要是包括rpc文件的遥感卫星影像，都可以利用rpc文件进行几何校正，尽管rpc几何校正有其局限之处。如果熟悉GDAL库，本文给出的代码很容易修改并应用到别的遥感卫星（例如国产高分一号、二号等卫星）。因此，本篇的博客内容具有较强的实际应用参考价值。

数据集

高分三号影像

关于高分三号 SAR卫星的介绍，网上有很多文章，这里不过多赘述了。

本文使用免费可用的数据高分3号数据(C波段全极化数据)进行介绍，免费可用的高分三号影像可以从PolSARpro官网样本数据进行下载。

本文使用上图中的San Francisco(美国旧金山)，Paris（法国巴黎）,Rennes（法国雷恩）三个地区的高分三号全极化SAR影像。

这三个地区中，旧金山地区临海，而雷恩市和巴黎市全位于陆地地区。

高分三号数据文件结构

简单介绍一下高分三号数据文件结构。下载的原始高分三号数据（全极化QPSI）解压得到数据文件集如下所示：

其中：

文件名	说明
GF3_KRN_QPSI_005777_E2.3_N48.8_20170914_L1A_AHV_L10002598957.incidence.xml	入射角数据文件
GF3_KRN_QPSI_005777_E2.3_N48.8_20170914_L1A_AHV_L10002598957.meta.xml	元数据文件
GF3_KRN_QPSI_005777_E2.3_N48.8_20170914_L1A_HH_L10002598957.jpg	HH通道jpg格式彩图
GF3_KRN_QPSI_005777_E2.3_N48.8_20170914_L1A_HH_L10002598957.rpc	HH通道对应的rpc文件
GF3_KRN_QPSI_005777_E2.3_N48.8_20170914_L1A_HH_L10002598957.thumb.jpg	HH通道缩略图jpg格式文件
GF3_KRN_QPSI_005777_E2.3_N48.8_20170914_L1A_HH_L10002598957.tiff	真正记录HH通道散射强度的tif文件(SLC)

DEM数据

DEM数据可以在cgiar网站https://srtm.csi.cgiar.org/srtmdata/下载得到，但需要科学上网。这里四个地区对应提供原始SRTMV4.1 DEM数据(分辨率为90m)，见下面的百度云盘链接：
链接：https://pan.baidu.com/s/1RPZ6I–hpnyfkLnuf9OvVw
提取码：63id

地区	DEM压缩文件名
San_Francisco（旧金山）	srtm_12_05.zip
Paris(巴黎)	srtm_37_03.zip
Rennes(雷恩)	srtm_36_03.zip

辐射定标处理

参考遥感事公众号的高分3号数据处理之PolSARpro文章的介绍。
首先，从遥感事公众号获取林科院的高分三号辐射定标程序，该程序需要ENVI的IDL中运行：

成功运行结果：

在GF-3压缩文件解压后的原始数据目录下会生成一个SLC文件夹：

SLC文件夹下会有ENVI .bin格式的二进制文件生成（可以ENVI、QGIS等软件打开）,这和PolSARpro 处理别的SAR数据生成的SLC文件基本是一样的：

RPC几何校正

基本原理

RPC模型，即有理函数(rational polynomial coefficient ,RPC)模型，是几何校正模型常用的一种模型，在通用的遥感教材中，通常都有所涉及，这里不过介绍。可以参考博客100天专业知识01—RPC系统校正

思路

辐射定标得到SLC目录数据文件为ENVI bin格式文件，需要读取GF-3原始数据集中的.rpc文件解析RPC域信息，解析其与GDAL库对应的RPC参数，并向.bin格式文件写入rpc域元数据，随后调用python中GDAL库内置的gdalwarp工具（即GDAL.Warp()函数）执行RPC校正，并保存为PolSARpro 对应ENVI格式的.bin文件（新的SLC目录数据集）。

实现

先看下.rpc文件的内容，以巴黎地区的数据集为例：

.rpc文字中的RPC参数与GDAL 官方定义的RPC参数基本是完全对应，仅仅是名字不同（这在使用程序解析提取需要注意的地方）。

在执行RPC校正请提前下载后对应的DEM数据
博主是在GF-3解压得到原数据目录下新建一个DEM文件夹，然后将对应的DEM zip文件放在该文件下，解压得到该DEM tif文件。

如果你按照我的文件夹组织方式处理的话，则后文中的python代码涉及路径需要改动的地方少一些。

陆地区域

需要提前安装好 python 的GDAL库，见python安装gdal的两种方法（最好使用第二种方法安装）。

使用python GDAL库实现的高分三号全极化SAR影像几何校正的代码如下（以巴黎地区为例）：

【版本信息：python 3.6 (python 3的版本应该都可以的)；GDAL V2.4.1（>1.8.0应该都是可以的）】

在运行代码之前，法国雷恩地区的数据文件集组织形式如下：

运行下面的代码后，将会解压后的原始数据目录下产生一个Geometric_Correction文件夹，其再下一级又有SLC和TEMP两个文件夹，其中SLC为几何校正后的S2复矩阵的ENVI bin文件，TEMP存放的是临时的bin格式文件（仅仅写入了RPC域元数据但没有校正的bin格式文件），如果觉得TEMP文件夹占用较多的存储空间，可以在RPC校正完成后将其删掉。

下面的代码沿着我的思路去看，加上我作的注释应该是比较好懂。不过，即使你看不懂，只需要注意路径问题和DEM问题，也可以很轻便地应用于你自身的GF-3数据集。

# -*- coding: utf-8 -*-
"""
-------------------------------------------------
# @Date     :2021/5/11 18:35
# @Author   :超级禾欠水
# @Email   : 1490930522@qq.com
-------------------------------------------------
"""import os
from osgeo import gdal
import glob"""
## 使用本代码请根据自身实际数据路径进行调整
## 有DEM数据情况，只需要修改origin_dir，old_SLC_dir，dem_tif_file 这三个路径参数
## 无DEM数据情况，只需要修改origin_dir, old_SLC_dir 这两个路径参数，注释掉dem_tif_file
"""# # Rennes,法国雷恩，全位于陆地
#origin_dir为高分三号原始数据解压的目录即.rpc文件目录
origin_dir = r'G:\test_GF3\2091067_Rennes'
# 获取原始数据目录下所有的rpc文件，实际上四个极化通道的rpc文件都是相同
rpc_files = list(glob.glob(os.path.join(origin_dir, '*.rpc')))# 使用林科院的GF-3辐射定标后得到的SLC目录
old_SLC_dir = r'G:\test_GF3\2091067_Rennes\SLC'# 存放几何校正后得到的文件目录
out_dir = os.path.join(origin_dir, 'Geometric_Correction')
if not os.path.exists(out_dir):os.mkdir(out_dir)# out_dir后创建一个存放中间转换的bin文件的TEMP临时目录，程序完成后可以删除整个TEMP临时目录
temp_bin_dir = os.path.join(out_dir, 'TEMP')
if not os.path.exists(temp_bin_dir):os.mkdir(temp_bin_dir)# out_dir后创建一个存放几何校正后的SLC下的S2复数矩阵的bin文件，以便后续在PolSARpro 中进行后续极化SAR处理
new_SLC_dir = os.path.join(out_dir, 'SLC')
if not os.path.exists(new_SLC_dir):os.mkdir(new_SLC_dir)# DEM tif文件的路径
# 注意DEM的覆盖范围要比原影像的范围大，此外，DEM不能有缺失值，有缺失值会报错
# 通常DEM在水域是没有值的（即缺失值的情况），因此需要将其填充设置为0，否则在RPC校正时会报错
# 这里使用的DEM是填充0值后的SRTM V4.1 3秒弧度的DEM(分辨率为90m)
# 请网上自行搜索下载GF-3 SAR影像对应的DEM文件（GeoTiff格式， WGS84坐标系）
dem_tif_file = r'G:\test_GF3\2091067_Rennes\DEM\srtm_36_03\srtm_36_03.tif'# 获取原始SLC目录下S2复数矩阵的的四个bin文件（s11.bin, s12.bin, s13.bin, s22.bin）
old_bin_files = list(glob.glob(os.path.join(old_SLC_dir, 's*.bin')))# ## San_Francisco，美国旧金山，有海域
# #origin_dir为高分三号原始数据解压的目录即.rpc文件目录，实际上四个极化通道的rpc文件都是相同
# origin_dir = r'G:\test_GF3\2599253_San_Francisco'
# # 获取原始数据目录下所有的rpc文件
# rpc_files = list(glob.glob(os.path.join(origin_dir, '*.rpc')))
#
# # 使用林科院的GF-3辐射定标后得到的SLC目录
# old_SLC_dir = r'G:\test_GF3\2599253_San_Francisco\SLC'
#
# # 存放几何校正后得到的文件目录
# out_dir = os.path.join(origin_dir, 'Geometric_Correction')
# if not os.path.exists(out_dir):
#     os.mkdir(out_dir)
#
# # out_dir后创建一个存放中间转换的bin文件的TEMP临时目录，程序完成后可以删除整个TEMP临时目录
# temp_bin_dir = os.path.join(out_dir, 'TEMP')
# if not os.path.exists(temp_bin_dir):
#     os.mkdir(temp_bin_dir)
#
# # out_dir后创建一个存放几何校正后的SLC下的S2复数矩阵的bin文件，以便后续在PolSARpro 中进行后续极化SAR处理
# new_SLC_dir = os.path.join(out_dir, 'SLC')
# if not os.path.exists(new_SLC_dir):
#     os.mkdir(new_SLC_dir)
#
# # DEM tif文件的路径
# # 注意DEM的覆盖范围要比原影像的范围大，此外，DEM不能有缺失值，有缺失值会报错
# # 通常DEM在水域是没有值的（即缺失值的情况），因此需要将其填充设置为0，否则在RPC校正时会报错
# # 这里使用的DEM是填充0值后的SRTM V4.1 3秒弧度的DEM(分辨率为90m)
# # 请网上自行搜索下载GF-3 SAR影像对应的DEM文件（GeoTiff格式， WGS84坐标系）
# dem_tif_file = r'G:\test_GF3\2599253_San_Francisco\DEM\srtm_12_05\srtm_12_05_fill_nodata.tif'
#
# # 获取原始SLC目录下S2复数矩阵的的四个bin文件（s11.bin, s12.bin, s13.bin, s22.bin）
# old_bin_files = list(glob.glob(os.path.join(old_SLC_dir, 's*.bin')))# 解析.rpc文件的函数
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# 将ENVI bin文件写入RPC域信息
def write_rpc_to_bin(bin_file, rpc_file):rpc_dict = parse_rpc_file(rpc_file)# 读取.bin文件old_ds = gdal.Open(bin_file)tif_X = old_ds.RasterXSizetif_Y = old_ds.RasterYSize# 获取像素数组array = old_ds.ReadAsArray()out_file = os.path.join(temp_bin_dir, os.path.basename(bin_file))# 创建临时的tif文件，并写入RPC域的文件envi_driver= gdal.GetDriverByName('ENVI')out_ds = envi_driver.CreateCopy(out_file, old_ds, 1)#out_ds.SetProjection('')# 向tif影像写入rpc域信息# 注意，这里虽然写入了RPC域信息，但实际影像还没有进行实际的RPC校正# 尽管有些RS/GIS能加载RPC域信息，并进行动态校正for k in rpc_dict.keys():out_ds.SetMetadataItem(k, rpc_dict[k], 'RPC')band= out_ds.GetRasterBand(1)band.WriteArray(array)out_ds.FlushCache()del old_ds, out_dsfor old_bin_file, rpc_file in zip(old_bin_files, rpc_files):write_rpc_to_bin(old_bin_file, rpc_file)print('成功向%s文件写入RPC域信息！'%(os.path.basename(old_bin_file)))# 提取刚写入rpc域信息的bin格式影像
uncorr_bins = glob.glob(os.path.join(temp_bin_dir, 's*.bin'))
# rpc校正后存放的tif文件，与old_bin_file对应，仅后缀名改为了.tif
corr_bins = [os.path.join(new_SLC_dir, os.path.basename(old_bin_file)) for old_bin_file in old_bin_files]## 设置rpc校正的参数
# 原图像和输出影像缺失值设置为0，输出影像坐标系为WGS84(EPSG:4326), 重采样方法为双线性插值（bilinear，还有最邻近‘near’、三次卷积‘cubic’等可选)
# RPC_DEM=G:\\test_GF3\\2599253_San_Francisco\\Geometric_Correction\\srtm_12_05_fill_nodata.tif 中
# G:\\test_GF3\\2599253_San_Francisco\\Geometric_Correction\\srtm_12_05_fill_nodata.tif为覆盖原影像范围的DEM
# 注意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代码
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"])## 对于全海域的影像或者不使用DEM校正的话，可以将transformerOptions有关的RPC DEM关键字删掉
## 即将上面gdal.WarpOptions注释掉，将下面的语句取消注释，无DEM时，影像范围的高程默认全为0
# wo = gdal.WarpOptions(srcNodata=0, dstNodata=0, dstSRS='epsg:4326',resampleAlg='bilinear', rpc=True, warpOptions=["INIT_DEST=NO_DATA"])# 耗费的时间与影像大小有关，请耐心等候
for corr_bin, uncorr_bin in zip(corr_bins, uncorr_bins):# 执行rpc校正wr = gdal.Warp(corr_bin,  uncorr_bin, options=wo)print('得到RPC校正后的%s影像'%os.path.basename(corr_bin))# 创建config.txt文件的函数
def create_config(dir, Nrow, Ncol):with open(os.path.join(dir, 'config.txt'), 'w') as f:f.write('Nrow\n')f.write('%d\n'%Nrow)f.write('--------\n')f.write('Ncol\n')f.write('%d\n'%Ncol)f.write('--------\n')f.write('PolarCase\n')f.write('bistatic\n')f.write('---------\n')f.write('PolarType\n')f.write('full\n')return# 创建config.txt文件
create_config(new_SLC_dir,wr.RasterYSize, wr.RasterXSize)del wr

运行效果（耗时多少与处理的数据集、电脑性能等有关，请耐心等候）：

生成新的S2矩阵目录（Geometric_Correction下的SLC）：

我不打算以这幅影像来检查校正效果，因为陆地相对而言对比不是很明显。可以初步看看几何效果，这是后面利用PolSARpro将其S2矩阵转换为T3矩阵（5x5多视）得到的PauliRGB图：

濒海区域

以美国的旧金山地区为例（srtm_12_05.tif），上面的代码需要修改的地方，主要是路径（将巴黎地区的路径代码注释掉，把旧金山地区的路径代码去掉注释），需要特别注意DEM文件：

注意：这里使用的DEM是用0值填充了缺失值（nodata）的DEM。你可以用没有填充缺失的DEM试下，程序报错，因为DEM的缺失值是-32768，这会导致RPC校正得到负的行列坐标值，从而报错而失败。

填充DEM缺失值

以旧金山地区对应的DEM zip文件为例：
（以下代码还需要numpy库的支持）

# -*- coding: utf-8 -*-
"""
-------------------------------------------------
# @Date     :2021/5/11 18:32
# @Author   :超级禾欠水
# @Email   : 1490930522@qq.com
-------------------------------------------------
""""""
## 填充临海地区的DEM影像缺失值
## 请结合自身实际研究区修改输入的DEM路径input_dem_tif
"""import os
from osgeo import gdal
import numpy as np# San_Francisco，旧金山 DEM
input_dem_tif = r'G:\test_GF3\2599253_San_Francisco\DEM\srtm_12_05\srtm_12_05.tif'output_dem_tif = os.path.join(os.path.dirname(input_dem_tif), os.path.basename(input_dem_tif)[:-4] + '_fill_nodata.tif')
old_ds = gdal.Open(input_dem_tif)
tif_X = old_ds.RasterXSize
tif_Y = old_ds.RasterYSize
array = old_ds.ReadAsArray()
geotrans = old_ds.GetGeoTransform()
proj = old_ds.GetProjection()
del old_ds# 将缺失值-32768填充为0
array = np.where(array==-32768, 0, array)# 将填充缺失值后的DEM保存为新的tif影像
gtiff_driver = gdal.GetDriverByName('GTiff')# 数据类型为有符号整型UInt16
out_ds = gtiff_driver.Create(output_dem_tif, tif_X, tif_Y, 1, gdal.GDT_UInt16)
out_ds.SetGeoTransform(geotrans)
out_ds.SetProjection(proj)
band = out_ds.GetRasterBand(1)
band.WriteArray(array)
out_ds.FlushCache()
del out_ds

运行成功后，会在原DEM 文件名下生成一个一个添加了后缀"_fill_nodata"的GeoTiff格式DEM 文件，可用于该区域后面的RPC校正：

RPC校正代码修改部分（即将雷恩地区部分的文件路径信息等代码注释，把美国旧金山地区部分的文件路径信息等代码取消注释）

# # # Rennes,法国雷恩，全位于陆地
# #origin_dir为高分三号原始数据解压的目录即.rpc文件目录
# origin_dir = r'G:\test_GF3\2091067_Rennes'
# # 获取原始数据目录下所有的rpc文件，实际上四个极化通道的rpc文件都是相同
# rpc_files = list(glob.glob(os.path.join(origin_dir, '*.rpc')))
#
# # 使用林科院的GF-3辐射定标后得到的SLC目录
# old_SLC_dir = r'G:\test_GF3\2091067_Rennes\SLC'
#
# # 存放几何校正后得到的文件目录
# out_dir = os.path.join(origin_dir, 'Geometric_Correction')
# if not os.path.exists(out_dir):
#     os.mkdir(out_dir)
#
# # out_dir后创建一个存放中间转换的bin文件的TEMP临时目录，程序完成后可以删除整个TEMP临时目录
# temp_bin_dir = os.path.join(out_dir, 'TEMP')
# if not os.path.exists(temp_bin_dir):
#     os.mkdir(temp_bin_dir)
#
# # out_dir后创建一个存放几何校正后的SLC下的S2复数矩阵的bin文件，以便后续在PolSARpro 中进行后续极化SAR处理
# new_SLC_dir = os.path.join(out_dir, 'SLC')
# if not os.path.exists(new_SLC_dir):
#     os.mkdir(new_SLC_dir)
#
# # DEM tif文件的路径
# # 注意DEM的覆盖范围要比原影像的范围大，此外，DEM不能有缺失值，有缺失值会报错
# # 通常DEM在水域是没有值的（即缺失值的情况），因此需要将其填充设置为0，否则在RPC校正时会报错
# # 这里使用的DEM是填充0值后的SRTM V4.1 3秒弧度的DEM(分辨率为90m)
# # 请网上自行搜索下载GF-3 SAR影像对应的DEM文件（GeoTiff格式， WGS84坐标系）
# dem_tif_file = r'G:\test_GF3\2091067_Rennes\DEM\srtm_36_03\srtm_36_03.tif'
#
# # 获取原始SLC目录下S2复数矩阵的的四个bin文件（s11.bin, s12.bin, s13.bin, s22.bin）
# old_bin_files = list(glob.glob(os.path.join(old_SLC_dir, 's*.bin')))
### San_Francisco，美国旧金山，有海域
#origin_dir为高分三号原始数据解压的目录即.rpc文件目录，实际上四个极化通道的rpc文件都是相同
origin_dir = r'G:\test_GF3\2599253_San_Francisco'
# 获取原始数据目录下所有的rpc文件
rpc_files = list(glob.glob(os.path.join(origin_dir, '*.rpc')))# 使用林科院的GF-3辐射定标后得到的SLC目录
old_SLC_dir = r'G:\test_GF3\2599253_San_Francisco\SLC'# 存放几何校正后得到的文件目录
out_dir = os.path.join(origin_dir, 'Geometric_Correction')
if not os.path.exists(out_dir):os.mkdir(out_dir)# out_dir后创建一个存放中间转换的bin文件的TEMP临时目录，程序完成后可以删除整个TEMP临时目录
temp_bin_dir = os.path.join(out_dir, 'TEMP')
if not os.path.exists(temp_bin_dir):os.mkdir(temp_bin_dir)# out_dir后创建一个存放几何校正后的SLC下的S2复数矩阵的bin文件，以便后续在PolSARpro 中进行后续极化SAR处理
new_SLC_dir = os.path.join(out_dir, 'SLC')
if not os.path.exists(new_SLC_dir):os.mkdir(new_SLC_dir)# DEM tif文件的路径
# 注意DEM的覆盖范围要比原影像的范围大，此外，DEM不能有缺失值，有缺失值会报错
# 通常DEM在水域是没有值的（即缺失值的情况），因此需要将其填充设置为0，否则在RPC校正时会报错
# 这里使用的DEM是填充0值后的SRTM V4.1 3秒弧度的DEM(分辨率为90m)
# 请网上自行搜索下载GF-3 SAR影像对应的DEM文件（GeoTiff格式， WGS84坐标系）
dem_tif_file = r'G:\test_GF3\2599253_San_Francisco\DEM\srtm_12_05\srtm_12_05_fill_nodata.tif'# 获取原始SLC目录下S2复数矩阵的的四个bin文件（s11.bin, s12.bin, s13.bin, s22.bin）
old_bin_files = list(glob.glob(os.path.join(old_SLC_dir, 's*.bin')))

运行结果：

生成的新的SLC目录：

注意事项

要用DEM文件进行RPC校正，DEM文件需满足的条件要求：

DEM文件坐标系为WGS84地理坐标系；
DEM文件范围覆盖整个遥感图像的范围，如果一个DEM文件没有覆盖完全，则需要镶嵌DEM文件（这里的数据集没有涉及到）；
DEM文件不能有缺失值；
此外，没有DEM的地区影像（如极地、海域等）也可以使用RPC来校正，但这时默认高程值为0。

后续处理

简单介绍一些RPC校正后得到GF-3 S2复矩阵文件在PolSARpro中的读取识别操作(以旧金山地区为例)。

S2矩阵

点击Save & Exit后，自动识别该SLC文件夹为S2矩阵目录环境，并打开GIMP产生mask_valid_pixels掩膜文件：

成功识别后，会有新的.bin.hdr等文件生成：

PolSARpro在坐标系上的bug

当你很兴奋地看到PolSARpro中成功识别S2矩阵元素的结果时，一个新的烦恼又跳出来。由于PolSARpro默认会生成新的.bin.hdr文件并且会使用默认坐标系，这将会使RPC校正得到.hdr头文件的坐标系信息在后续操作时并忽略。

例如，当你用文本编辑器打开s11.hdr和s11.bin.hdr时，你会发现两者的内容中，
map_info信息是不一样的。（注意，coordinate system string信息对于ENVI bin文件不是必要的，但map_info的信息不应缺少）。

修正这一坐标的方式，就是将RPC校正得到.hdr文件中正确的map_info信息替换掉PolSARpro默认生成的错误的.bin.hdr的map_info信息。考虑.bin.hdr文件数量较多，
这里使用python写了一个简短的更正.bin.hdr头文件坐标信息代码，完成该操作。
（结合自身实际情况修改代码中的路径）

# -*- coding: utf-8 -*-
"""
-------------------------------------------------
# @Date     :2021/5/13 11:38
# @Author   :超级禾欠水
# @Email   : 1490930522@qq.com
-------------------------------------------------
"""import os
import glob# 旧金山地区RPC校正后得到的SLC目录
new_SLC_dir = r'G:\test_GF3\2599253_San_Francisco\Geometric_Correction\SLC'# # 法国雷恩地区RPC校正后得到的SLC目录
# new_SLC_dir = r'G:\test_GF3\2091067_Rennes\Geometric_Correction\SLC'true_hdr_file = os.path.join(new_SLC_dir, 's11.hdr')
with open(true_hdr_file) as f:text = f.readlines()
# 提取坐标系统信息
map_info = text[-4]# 要修改含有错误坐标系统的.bin.hdr文件的目录路径
bin_hdr_dir = r'G:\test_GF3\2599253_San_Francisco\Geometric_Correction\SLC'
error_bin_hdrs = glob.glob(os.path.join(bin_hdr_dir, '*.bin.hdr'))
for error_bin_hdr in error_bin_hdrs:with open(error_bin_hdr, 'r') as fi:lines =fi.readlines()for i, line in enumerate(lines):if  'map info' in line:lines[i] = map_infowith open(error_bin_hdr, 'w') as fo:fo.writelines(lines)
print('成功修改.bin.hdr文件的坐标信息!')

运行结果：

可以看到.bin.hdr的map_info信息与.hdr文件的对应一致了：

这个bug不仅发生在识别S2矩阵时，在生成T3等矩阵元素也可能发生。若后者情况发生，也同样可以替换掉map_info信息的办法进行处理。

校正效果检验

以旧金山地区为例，因为我对这个区域相对较熟悉一些。
未几何校正前的快视图(对应的.thumb.jpg后缀文件)，这里是属于上下颠倒的情况，说明这时卫星是右侧视升轨状态：

这里因为Google Earth在国内被禁了，这里使用旧金山的Sentinel-2光学影像来检查GF-3的几何校正效果。GF-3的分辨率为8m,Sentinel-2 的分辨率为10m。

如果你自身的高分三号影像覆盖区域与本文不一样，与这里相似，可以自己下载对应区域相近时间段的Sentinel-2等光学卫星影像来检查其RPC校正效果。

如下图所示：

S11.bin为已进行RPC校正后得到GF-3 S2 复矩阵S11元素（HH通道）影像（2017年9月15日）,分辨率为8m；

S2A_20170917_San_Francisco_WGS84.tif为旧金山Sentinel-2A 真彩色合成影像（2017年9月24日）, 分辨率为24m, 已转为WGS84坐标系。尽管这幅影像有点云，但不影响检查校正效果。

原Sentinel-2A文件名为
S2A_MSIL1C_20170917T190351_N0205_R113_T10SEG_20170917T190345.SAFE，Sentinel卫星影像的下载不必再说了。文件名S2A_20170917_San_Francisco_WGS84.tif百度云盘下载如下：

链接：https://pan.baidu.com/s/17Sw5VyqwAw98drrToe_APg
提取码：ehvr

借助ENVI工具来检查一下（直接将RPC校正后得到的s12.bin.hdr或是s12.hdr拖入ENVI中即可成功读取）：

注意，观测海陆边界，可以将看到两者是基本吻合，证实了本文代码实现的RPC几何校正的有效性和可靠性。上方框海陆边界吻合得很好，下方框海陆边界则有点区别，这主要是SAR和光学卫星成像方式不一样导致的。（放大后可以粗略判断RPC校正后得到的GF-3影像与Sentinel-2影像的配准误差不超过2个像素，对于一般的分类等应用问题不大。如果对精度要求更高的话，可以在RPC校正结果的基础上进一步进行配准。）

T3矩阵

接上面的操作，按下图设置转换参数，点击Run即可：

成功后的界面，会弹出有PauliRGB.bmp或mask_valid_pixels.bmp：

得到了如上图的非常漂亮的PauliRGB图，并有新生成的T3文件夹：

后续的利用PolSAR处理GF-3 全极化数据的教程请参考遥感事公众号高分3号数据处理之PolSARpro等文章
或请参考我之前写的利用PolSARpro处理Sentinel-1卫星SAR影像的几篇博客。

另外1个地区GF-3影像不做演示，留给各位测试练习。

后语

本人后面的研究方向不在极化SAR影像应用方面了，后面创作的博客将很少讨论极化SAR的内容。其实可以做极化SAR处理的工具不只PolSARpro一个开源工具，还有OTB(Orfeo Toolbox)等开源工具，有兴趣的同志可以研究一下，后面可能会介绍OTB工具的使用。后面精力主要放在科研上（论文），更新博客的速度后慢一些，不过，我努力写出一点有价值的内容！谢谢关注，祝好！

参考文献

[1] 博园客博客文章使用GDAL工具对OrbView-3数据进行正射校正
[2] GDAL库gdalwarp工具document. https://gdal.org/programs/gdalwarp.html
[3] 李民录. GDAL源码剖析与开发指南[M]. 人民邮电出版社, 2014.
[4] 李鹏, 黎达辉, 李振洪*, 王厚杰. 黄河三角洲地区GF-3雷达数据与Sentinel-2多光谱数据湿地协同分类研究[J]. 武汉大学学报(信息科学版), 2019, 44(11):1641-1649.

使用Python GDAL库对高分三号全极化SAR影像进行RPC几何校正（PolSARpro格式）相关推荐

python 读取geotiff_科学网—利用python GDAL库读写geotiff格式的遥感影像方法 - 张伟的博文...
(1)利用python GDAL库读写geotiff格式的遥感影像方法,具有很好的参考价值,不错! from osgeo import gdal import numpy as np def read ...
python读取tiff影像_科学网—利用python GDAL库读写geotiff格式的遥感影像方法 - 张伟的博文...
(1)利用python GDAL库读写geotiff格式的遥感影像方法,具有很好的参考价值,不错! from osgeo import gdal import numpy as np def read ...
高分三号（GF-3）单极化/双极化/全极化多时相变化检测
高分三号SAR成像卫星具备聚束单极化(1m).超精细条带单极化(3m).精细条带1双极化(5m)和全极化条带1(8m)等12种成像模式. 雷达影像信息处理系统SARInfo具备高分三号单极化.双极化. ...
[转载] python标准库系列教程(三)——operator库详细教程
参考链接: Python中的Inplace运算符| 2(ixor(),iand(),ipow()等) python进阶教程机器学习深度学习长按二维码关注进入正文 Python基础学习:oper ...
1000道python题库系列_1000道Python题库系列分享三
上一期题目链接:1000道python题库系列分享二(48道) 上一题题目参考答案: 2.1 31 2.2 'F' 2.3 python采用的是基于值得内存管理方式,在Python中可以为不同变量赋值 ...
高分二号多光谱GEOTIFF影像批量转换为JPG图片
最近在做语义分割方面的研究,想对比一下高分二号4个波段的模型输入和JPG输入的结果有什么差别.使用了一些图片转换软件发现转换后的jpg图片一片漆黑,应该是geitif格式的特殊性吧!所有在网上找了一些 ...
Python gdal库读取tif文件
from osgeo import gdal # GDAL库主要提供对栅格数据的处理,使用抽象数据模型来解析所支持的数据格式 import filename_cut as fc import matp ...
【python】Romantic Rose(玫瑰三号)
# 第三种:画法的玫瑰没有第二种的好看 import turtledef init():turtle.setup(width=0.9, height=0.9)turtle.speed(10)def f ...
基于GDAL库图像读写——涉及(tif/tiff/bmp/jpg/png/gif等)多种格式图像的I/O
说来也是神奇,之所以会发这篇文章,还是因为在CSDN发文章中无法上传tif格式的图片造成的,如下. 转入文件夹下,并不显示tif图片,选择所有文件确认一副tiff图片,仍然提示上传错误.不知道CSDN ...
Python扩展库numpy中where()函数的三种用法
第一种用法:只给where()函数传递一个数组作为参数,返回其中非0元素的下标. 第二种用法:给where()函数传递一个包含True/False值的数组,返回该数组中True值的下标,结合numpy ...

使用Python GDAL库对高分三号全极化SAR影像进行RPC几何校正（PolSARpro格式）

使用Python GDAL库对高分三号全极化SAR影像进行RPC几何校正（PolSARpro格式）

前言

数据集

高分三号影像

高分三号数据文件结构

DEM数据

辐射定标处理

RPC几何校正

基本原理

思路

实现

陆地区域

濒海区域

填充DEM缺失值

注意事项

后续处理

S2矩阵

PolSARpro在坐标系上的bug

校正效果检验

T3矩阵

后语

参考文献

使用Python GDAL库对高分三号全极化SAR影像进行RPC几何校正（PolSARpro格式）相关推荐

最新文章

热门文章