ICode9
精准搜索请尝试： 精确搜索
首页
编程语言>
数据库
系统相关
互联网
其他分享
Python
Java
JavaScript
android
PHP
首页 > 其他分享> 文章详细

FY-4A批量几何校正（基于GDAL）

2021-06-10 22:05:21  阅读：489  来源： 互联网

标签：校正 FY gdal EPSG -_ 4A tif GDAL Tablefile


引言

​ 前次介绍了FY-4A遥感影像的几何校正方法，然而其主体基于ENVI实现，虽然可以通过ENVI_DO_IT和IDL语言实现批量化处理，但由于ENVI中GLT的生成需要经纬度查找表没有无效值，需要裁剪，所以比较麻烦。

​ 下面介绍一下基于gdal进行FY-4A几何校正的方法。仔细阅读李民录老师《gdal实现静止卫星几何校正》相关文章后，发现其主要流程都是先gdalinfo查看数据集，再gdaltranslate生成vrt，然后修改vrt添加查找表数据路径，最后gdalwarp进行几何校正。

​ 实施过程中也发现了一些问题，如圆盘数据的类型不同，有HDF/NC/Tiff，存储方式不同，读取时文件的路径也各不相同，有的文件类型如NC甚至无法生成VRT文件。综合分析后，笔者决定拆解这一问题，既然gdalwarp的核心是根据影像数组，经纬度查找表数组进行几何校正，那就可以只凭借这三个数组实现几何校正，故可以生产各自的geotiff文件，略去多余的信息。

步骤

  1. 数据准备

    1.1 遥感影像nc转geotiff

    ​ FY-4A遥感影像属于netCDF数据集，在使用gdal_translate.exe生成vrt文件时，netCDF文件由于库自身的问题会导致错误，故首先需要把FY-4A数据转换成GeoTiff数据。

​ 1.2 生成lon/lat对应的geotiff数据

​ 即经纬度查找表

  1. gdal_translate生成VRT文件

    有关gdal的命令如gadl_translate.exe/gdalwarp.exe可以在cmd里调用，也可以在python环境下调用gdal函数。（需要注意在cmd里调用需要声明插件的路径（proj.db/...），或者在环境变量中声明。）

gdal_translate.exe -of VRT E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif E:/Tablefile/before.vrt

gdal.Translate('E:/Tablefile/before.vrt', 'E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif', options=gdal.TranslateOptions(format='VRT'))

  1. 修改VRT文件，添加用于校正的lon/lat数据集

    以上操作生成VRT文件(beforer.vrt)内容如下：

<VRTDataset rasterXSize="2748" rasterYSize="2748">
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>-999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2">
    <NoDataValue>-999</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>2</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="3">
    <NoDataValue>-999</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>3</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

遥感影像有几个波段，就会生成几个VRTRasterBand，部分数据（单波段、部分多波段）不会生成"TIFFTAG_XRESOLUTION"信息（情况暂时不明），需要自行添加，重复添加不会有影响。（记得在SourceFilename处添加波段路径）

本次校正修改如下：

添加经纬度查找表数据描述信息（在X_DATASET/Y_DATASET处添加lon/lat文件路径）

Modis数据由于自带lon/lat数据集，所以会自动生成这段信息；而其他不含lon/lat的卫星数据（如FY-4A、单独的geotiff）需要手动添加这段信息。

  <Metadata>
    <MDI key="TIFFTAG_XRESOLUTION">1</MDI>
    <MDI key="TIFFTAG_YRESOLUTION">1</MDI>
  </Metadata>
  <Metadata domain="GEOLOCATION">
    <MDI key="LINE_OFFSET">1</MDI>
    <MDI key="LINE_STEP">1</MDI>
    <MDI key="PIXEL_OFFSET">1</MDI>
    <MDI key="PIXEL_STEP">1</MDI>
    <MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]</MDI>
    <MDI key="X_BAND">1</MDI>
    <MDI key="X_DATASET">E:\Tablefile\CBHRegression\geometryCorrect\lookuptable\4km\lonDisk.tif</MDI>
    <MDI key="Y_BAND">1</MDI>
    <MDI key="Y_DATASET">E:\Tablefile\CBHRegression\geometryCorrect\lookuptable\4km\latDisk.tif</MDI>
  </Metadata>

修改后VRT文件(after.vrt)如下：

此处可以更改dataType，设置数据类型。

<VRTDataset rasterXSize="2748" rasterYSize="2748">
  <Metadata>
    <MDI key="TIFFTAG_XRESOLUTION">1</MDI>
    <MDI key="TIFFTAG_YRESOLUTION">1</MDI>
  </Metadata>
  <Metadata domain="GEOLOCATION">
    <MDI key="LINE_OFFSET">1</MDI>
    <MDI key="LINE_STEP">1</MDI>
    <MDI key="PIXEL_OFFSET">1</MDI>
    <MDI key="PIXEL_STEP">1</MDI>
    <MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]</MDI>
    <MDI key="X_BAND">1</MDI>
    <MDI key="X_DATASET">E:\Tablefile\CBHRegression\geometryCorrect\lookuptable\4km\lonDisk.tif</MDI>
    <MDI key="Y_BAND">1</MDI>
    <MDI key="Y_DATASET">E:\Tablefile\CBHRegression\geometryCorrect\lookuptable\4km\latDisk.tif</MDI>
  </Metadata>
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>-999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2">
    <NoDataValue>-999</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>2</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="3">
    <NoDataValue>-999</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>3</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>
  1. gdalwarp几何校正

最后利用gdalwarp进行几何校正即可（推荐使用函数进行）

注意指定geoloc利用查找表进行校正、利用t_srs添加投影信息，这样会自动生成含有geotrans/proj的geotiff。

gdalwarp.exe -geoloc -t_srs EPSG:4326 E:/Tablefile/after.vrt E:/Tablefile/warp.tif
# 但是这种方式有时会找不到插件路径，需要自行设置

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
gdal.Warp(r'E:/Tablefile/warp.tif',
                r'E:/Tablefile/after.vrt',
                format='GTiff', geoloc=True,
                dstSRS=srs.ExportToWkt(), resampleAlg=gdal.GRIORA_NearestNeighbour)

结果

几何校正结果如下，黄色是-180180，-9090的全球地图框线，红色是东海附近框线。由于FY-4A只拍到了地球的一侧，所以只有右半区有成像，左半边的黑色区域对应圆盘的最右侧，其经度起于最左侧地区。（此时中国正在夜间）

后记

  1. 几何校正的质量和经纬度查找表的质量直接相关，错误的经纬度查找表会导致校正失败。（注意查找表的数据格式，必须是单波段的float32数组）

  2. FY-4A有HDF，也有NC数据集，转换TIFF时需要注意；此外部分数据如CTH读取后数组需要.astype(np.float32)，否则在gdal写入栅格数据时gdal内部会发生数据类型转换导致错误。

  3. 中国区REGC与全圆盘DISK操作方法相同，注意选择合适的经纬度查找表即可。

  4. 该方法可以在已知经纬度查找表数据、遥感影像数据的情况下进行几何校正。将以上过程代码化，即可完成遥感影像的批量几何校正。

标签：校正,FY,gdal,EPSG,-_,4A,tif,GDAL,Tablefile
来源： https://www.cnblogs.com/ljwgis/p/14872976.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

关于我们 | 联系我们 | 留言反馈

专注分享技术，共同学习，共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有