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

水深无声 2023-01-20 06:53 1468阅读 0赞

### 使用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官网[样本数据][Link 1]进行下载。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70]  
本文使用上图中的San Francisco(美国旧金山)，Paris（法国巴黎）,Rennes（法国雷恩）三个地区的高分三号全极化SAR影像。

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

### 高分三号数据文件结构 ###

简单介绍一下高分三号数据文件结构。下载的原始高分三号数据（全极化QPSI）解压得到数据文件集如下所示：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 1]  
其中：

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

## DEM数据 ##

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

<table> 
 <thead> 
  <tr> 
   <th>地区</th> 
   <th>DEM压缩文件名</th> 
  </tr> 
 </thead> 
 <tbody> 
  <tr> 
   <td>San_Francisco（旧金山）</td> 
   <td>srtm_12_05.zip</td> 
  </tr> 
  <tr> 
   <td>Paris(巴黎)</td> 
   <td>srtm_37_03.zip</td> 
  </tr> 
  <tr> 
   <td>Rennes(雷恩)</td> 
   <td>srtm_36_03.zip</td> 
  </tr> 
 </tbody> 
</table>

# 辐射定标处理 #

参考遥感事公众号的[高分3号数据处理之PolSARpro][3_PolSARpro]文章的介绍。  
首先，从遥感事公众号获取林科院的高分三号辐射定标程序，该程序需要ENVI的IDL中运行：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 2]  
成功运行结果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 3]  
在GF-3压缩文件解压后的原始数据目录下会生成一个SLC文件夹：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 4]  
SLC文件夹下会有ENVI .bin格式的二进制文件生成（可以ENVI、QGIS等软件打开）,这和PolSARpro 处理别的SAR数据生成的SLC文件基本是一样的：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 5]

# RPC几何校正 #

## 基本原理 ##

RPC模型，即有理函数(rational polynomial coefficient ,RPC)模型，是几何校正模型常用的一种模型，在通用的遥感教材中，通常都有所涉及，这里不过介绍。可以参考博客[100天专业知识01—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文件的内容，以巴黎地区的数据集为例：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 6]  
.rpc文字中的RPC参数与GDAL 官方定义的[RPC参数][RPC]基本是完全对应，仅仅是名字不同（这在使用程序解析提取需要注意的地方）。

**在执行RPC校正请提前下载后对应的DEM数据**  
博主是在GF-3解压得到原数据目录下新建一个DEM文件夹，然后将对应的DEM zip文件放在该文件下，解压得到该DEM tif文件。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 7]  
如果你按照我的文件夹组织方式处理的话，则后文中的python代码涉及路径需要改动的地方少一些。

### 陆地区域 ###

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

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

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

在运行代码之前，法国雷恩地区的数据文件集组织形式如下：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 8]

*运行下面的代码后，将会解压后的原始数据目录下产生一个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

运行效果（耗时多少与处理的数据集、电脑性能等有关，请耐心等候）：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 9]  
生成新的S2矩阵目录（Geometric\_Correction下的SLC）：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 10]

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

### 濒海区域 ###

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

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 12]

#### 填充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校正：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 13]

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')))

运行结果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 14]

生成的新的SLC目录：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 15]

## 注意事项 ##

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

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

# 后续处理 #

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

## S2矩阵 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 16]

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 17]

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 18]

### PolSARpro在坐标系上的bug ###

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

例如，当你用文本编辑器打开s11.hdr和s11.bin.hdr时，你会发现两者的内容中，  
map\_info信息是不一样的。（注意，coordinate system string信息对于ENVI bin文件不是必要的，但map\_info的信息不应缺少）。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 19]

修正这一坐标的方式，就是将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文件的坐标信息!')

运行结果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 20]

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 21]

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

### 校正效果检验 ###

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

这里因为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][https_pan.baidu.com_s_17Sw5VyqwAw98drrToe_APg]  
提取码：ehvr

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

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

## T3矩阵 ##

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 24]

成功后的界面，会弹出有PauliRGB.bmp或mask\_valid\_pixels.bmp：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 25]

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 26]

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

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

# 后语 #

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

# 参考文献 #

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

[Link 1]: https://ietr-lab.univ-rennes1.fr/polsarpro-bio/gaofen-3/
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70]: /images/20221021/8b85e2ad836b4a77aac8882a0bf3a477.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 1]: /images/20221021/2c953db564bc4864bed50d75084ec9ce.png
[https_srtm.csi.cgiar.org_srtmdata]: https://srtm.csi.cgiar.org/srtmdata/
[https_pan.baidu.com_s_1RPZ6I_hpnyfkLnuf9OvVw]: https://pan.baidu.com/s/1RPZ6I--hpnyfkLnuf9OvVw
[3_PolSARpro]: https://mp.weixin.qq.com/s?__biz=MzI0MDQ2NTU2Ng==&mid=2247483656&idx=1&sn=868780ababdbc69f2e72ccc4d8364302&chksm=e91b22f9de6cabefd926befddca6094781ab5a41231f90da19141b0252fbb561dcec5f6e9164&scene=126&sessionid=1620731897&key=65b4b4f73f1b09ad5e223739dfc6f54aa0b5c32dbb60c4fce7f8788a664b2dca0fae01779b619c6aa081076d7cf69d25a7b06735b4918d58815ac9a6c516584f137a9a1ad52a729c453f40bc8a757b653220acd907b2a62b69a87da03160cca3872d365b874f869f42fc3426f129a97dc3f402af312b076b2b65619c9b9fc31b&ascene=1&uin=MTE4OTUxNjIxNw==&devicetype=Windows%2010%20x64&version=6300002f&lang=zh_CN&exportkey=A4ZXkbSnBTya1DEtedBI130=&pass_ticket=HuxC4fqOK6laNS7Dr2lJJumB63DAD1cbJFRP9IYkMKPpNJzuCpgWGwUOTlPzCpJv&wx_header=0
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 2]: /images/20221021/1e368fb0f47f4846bfafe0ef9713a9b9.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 3]: /images/20221021/65ffdc6378c649fea219391ce3d4038b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 4]: /images/20221021/2c10dfe4bebc48db8fd2bdb9f4534981.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 5]: /images/20221021/908c43d024584d31bd3a6258b12958ca.png
[100_01_RPC]: http://blog.sina.com.cn/s/blog_ae72410c0102vl98.html
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 6]: /images/20221021/294383cb5025486ca30cc3da1f30aa56.png
[RPC]: http://geotiff.maptools.org/rpc_prop.html
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 7]: /images/20221021/de23c01b4d08464f88e0d046c8dc1b72.png
[python_gdal]: https://www.jb51.net/article/173083.htm
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 8]: /images/20221021/19f054f65aef4fd587820aef1f998b96.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 9]: /images/20221021/bff248e2ecab497d99976a91b69667a6.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 10]: /images/20221021/acb4615058564e55a48b143351d52937.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 11]: /images/20221021/65eb5b182ee94f859c27708f862eaf88.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 12]: /images/20221021/866ffc80241942f1a4da7b7df5e91a70.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 13]: /images/20221021/e79e2fb164844f72b071b6b30a61aae2.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 14]: /images/20221021/87bef79d7c024325858d7747786a49b5.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 15]: /images/20221021/97b7e62e79224ce9b14e2797c9b9d99c.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 16]: /images/20221021/33e63f55bbb54bbf86c27f3f1025c26e.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 17]: /images/20221021/86050b11dec144f7a21977132e4712ec.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 18]: /images/20221021/9c161d42435c4a56aaa37f28ab05a05b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 19]: /images/20221021/0ea9a28ba88e4e159b6756bf4630cb23.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 20]: /images/20221021/660960c080234dd9886ee7c91d11b0c7.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 21]: /images/20221021/366d91291a434ba987c9a81d72e29e92.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 22]: /images/20221021/13661d452b124234bbbbbcd66a159318.png
[https_pan.baidu.com_s_17Sw5VyqwAw98drrToe_APg]: https://pan.baidu.com/s/17Sw5VyqwAw98drrToe_APg
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 23]: /images/20221021/7cacb36ce49c4bfe992e19e49988f942.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 24]: /images/20221021/ce14b8d2e716462692f1f752663ec118.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 25]: /images/20221021/917eafa4ff81472e9ea94691b043b630.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xpZGFodWlsaWRhaHVp_size_16_color_FFFFFF_t_70 26]: /images/20221021/79270257bb5543aa87180c5c3bc9e64e.png
[3_PolSARpro 1]: https://mp.weixin.qq.com/s?__biz=MzI0MDQ2NTU2Ng==&mid=2247483656&idx=1&sn=868780ababdbc69f2e72ccc4d8364302&chksm=e91b22f9de6cabefd926befddca6094781ab5a41231f90da19141b0252fbb561dcec5f6e9164&scene=126&sessionid=1620804258&key=a614a7a47c2e2fdbcacdbbe8cc4b79d16b8ed18fa4079a946fbc02e2508d98f7239873f814b63ae48d9c9dc9aadeee13be56489859470a21bf410aeb4c0cd56e60af2bded0c03acda33313fcb51e069d8e934c1fac50024e713adce9486e4eccac67d6952ab8e7b7e996ce7173eea0ca65eb2726c9d30c710c22f1e7938bfbcd&ascene=1&uin=MTE4OTUxNjIxNw==&devicetype=Windows%2010%20x64&version=6300002f&lang=zh_CN&exportkey=AxSkhDxH9PBOnZB1MlsE6Go=&pass_ticket=NDwMJTdWKPL6qqZ3T1tomPzhdhoqCKE9X5iFIpvZIN43JVaEm73RaNKoGOK1MOwh&wx_header=0
[GDAL_OrbView-3]: https://www.cnblogs.com/xiaowangba/archive/2012/11/17/6314000.html
[https_gdal.org_programs_gdalwarp.html]: https://gdal.org/programs/gdalwarp.html