使用Python GDAL库对高分三号全极化SAR影像进行RPC几何校正（PolSARpro格式）
（原创文章，码字不易，转载请注明来源，谢谢！）
前言
很久没有更新博客，整个研三来基本都在忙着各种事情（考博、写毕业论文、找工作、党建等），现在基本收尾了。在毕业答辩之前，终于零星有点空余时间写下这篇博客，希望对国产卫星的应用推广有所帮助。
尽管高分三号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.html
rpc_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] = value
for 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.RasterXSize
tif_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_ds
for 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_info
with 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.
7524
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
