FY-4A批量幾何校正（基於GDAL）

本文轉載自查看原文 2021-06-10 21:47 1593 遙感數據介紹/ 遙感圖像處理算法

引言

前次介紹了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 遙感影像nc轉geotiff

FY-4A遙感影像屬於netCDF數據集，在使用gdal_translate.exe生成vrt文件時，netCDF文件由於庫自身的問題會導致錯誤，故首先需要把FY-4A數據轉換成GeoTiff數據。

1.2 生成lon/lat對應的geotiff數據

即經緯度查找表

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

修改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>

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)

結果

幾何校正結果如下，黃色是-180_180，-9090的全球地圖框線，紅色是東海附近框線。由於FY-4A只拍到了地球的一側，所以只有右半區有成像，左半邊的黑色區域對應圓盤的最右側，其經度起於最左側地區。（此時中國正在夜間）

后記

幾何校正的質量和經緯度查找表的質量直接相關，錯誤的經緯度查找表會導致校正失敗。（注意查找表的數據格式，必須是單波段的float32數組）
FY-4A有HDF，也有NC數據集，轉換TIFF時需要注意；此外部分數據如CTH讀取后數組需要.astype(np.float32)，否則在gdal寫入柵格數據時gdal內部會發生數據類型轉換導致錯誤。
中國區REGC與全圓盤DISK操作方法相同，注意選擇合適的經緯度查找表即可。
該方法可以在已知經緯度查找表數據、遙感影像數據的情況下進行幾何校正。將以上過程代碼化，即可完成遙感影像的批量幾何校正。

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 Python 利用GDAL對圖像進行幾何校正 PIE-Basic 幾何校正 ENVI【遙感圖像預處理之圖像的幾何校正】 ENVI軟件基本操作——輻射校正、輻射定標、大氣校正、幾何校正利用多項式實現圖像幾何校正(Matlab實現) 【OpenCV】幾何失真校正與局部增強對極幾何 GDAL------安裝GDAL 顏色校正 gamma校正