当前位置：首页 > 代码相关 > 正文内容

GDAL使用（五）

admin4年前 (2021-02-14)代码相关4370

在之前说的GDALWarp.exe程序的参数中有个校正选项，\[-order n | -tps | -rpc | -geoloc\]，这里说一下geoloc校正。

这里的内容还是参考李民录老师的书以及GDAL官方手册内容。

geoloc校正类似于ENVI的GLT校正，ENVI中就是建立一个地理查找表，通过地理查找表将没有任何空间参考的原始影像变换到具有合适坐标系的影像。

可以使用geoloc校正的数据有MODIS、风云卫星（FY）、海洋卫星（HY）等，我只用过MODIS和HY-1C数据，这里也以他们为例子做个记录。

这里的MODIS数据是MOD021KM产品，在powershell中使用gdalinfo来查看信息：

.\gdalinfo.exe E:\testdataMOD021KM.A2016352.0250.061.2017329005547.hdf

![modisinfo.jpg](https://lovedm.club/zb_users/upload/2021/02/202102141613314270892742.jpg "modisinfo.jpg")

可以在输出的子数据集中看到，图中红框内是经纬度波段。MODIS经纬度数据和图像数据不是一一对应的，所以大小是比原始影像大小要小的。

选取EV\_1KM\_RefSB数据集进行校正。在powershell中输入以下内容：

.\gdalwarp.exe -geoloc -t\_srs \'+proj=longlat +datum=WGS84 +no\_defs\' HDF4\_EOS:EOS\_SWATH:\"E:\testdataMOD021KM.A2016352.0250.061.2017329005547.hdf\":MODIS\_SWATH\_Type\_L1B:EV\_1KM\_RefSB E:\GEOREFMODIS.tif

工具会完成图像的校正，这里参数的意思是指定使用geoloc方式校正，输出坐标系采用WGS84，输入文件的路径为数据在数据集中的路径。最后是输出文件路径。

当然也可以直接校正到投影坐标系下，如下面：

.\\gdalwarp.exe -geoloc -t\_srs \'+proj=aea +lat\_1=15 +lat\_2=65 +lat\_0=30 +lon\_0=95 +x\_0=0 +y\_0=0 +datum=WGS84 +units=m +no\_defs\' HDF4\_EOS:EOS\_SWATH:"E:\testdata\MOD021KM.A2016352.0250.061.2017329005547.hdf":MODIS\_SWATH\_Type\_L1B:EV\_1KM\_RefSB E:\GEOREFMODIS2.tif

MODIS是以hdf4格式存储的，且GDAL工具能够识别出经纬度信息建立起查找表，所以直接使用工具进行校正即可。也可以和上篇文章一样在代码中调用gdalwarp函数进行校正，无非是换一下参数。

海洋卫星数据，这里以海洋一号C卫星的COCTS传感器L1B产品数据为例，分发下来的数据格式为hdf5，文件格式后缀以.h5结尾，这个是无法校正的，按照上面的方法校正时会出错，如下：

![COCTS校正error.jpg](https://lovedm.club/zb_users/upload/2021/02/202102141613314299450267.jpg "COCTS校正error.jpg")

错误大概是无法建立地理查找数组，那么怎么办呢，这里可以使用一个中间格式进行一下转换，这个格式是.vrt格式，这是一种虚拟格式，可以使用xml来进行存储，GDAL还有相应的工具构建vrt文件，当然，这里自己手动进行vrt文件的创建。

首先可以使用gdalinfo工具来查看数据集的信息，如下：

![COCTSinfo.jpg](https://lovedm.club/zb_users/upload/2021/02/202102141613314319705360.jpg "COCTSinfo.jpg")

从上面可以得到的可用信息为：

数据路径：以红光波段为例：HDF5:"E:\testdataCOCTSOrigin.h5"://Geophysical\_Data/L\_670

纬度路径：HDF5:"E:\testdata\COCTSOrigin.h5"://Navigation\_Data/Latitude

经度路径：HDF5:"E:\testdata\COCTSOrigin.h5"://Navigation\_Data/Longitude

然后构造vrt文件：

``` markup
<VRTDataset rasterXSize="1656" rasterYSize="2016">
  <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">HDF5:"E:/testdata/COCTSOrigin.h5"://Navigation_Data/Longitude</MDI>
   <MDI key="Y_BAND">1</MDI>
   <MDI key="Y_DATASET">HDF5:"E:/testdata/COCTSOrigin.h5"://Navigation_Data/Latitude</MDI>
  </Metadata>

<VRTRasterBand dataType="Float32" band="1">
     <SimpleSource>
       <SourceFilename relativeToVRT="1">HDF5:"COCTSOrigin.h5"://Geophysical_Data/L_412</SourceFilename>
       <SourceBand>1</SourceBand>
     </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2">
     <SimpleSource>
       <SourceFilename relativeToVRT="1">HDF5:"COCTSOrigin.h5"://Geophysical_Data/L_443</SourceFilename>
       <SourceBand>1</SourceBand>
     </SimpleSource>
  </VRTRasterBand> 
  <VRTRasterBand dataType="Float32" band="3">
     <SimpleSource>
       <SourceFilename relativeToVRT="1">HDF5:"COCTSOrigin.h5"://Geophysical_Data/L_490</SourceFilename>
       <SourceBand>1</SourceBand>
     </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="4">
     <SimpleSource>
       <SourceFilename relativeToVRT="1">HDF5:"COCTSOrigin.h5"://Geophysical_Data/L_520</SourceFilename>
       <SourceBand>1</SourceBand>
     </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="5">
     <SimpleSource>
       <SourceFilename relativeToVRT="1">HDF5:"COCTSOrigin.h5"://Geophysical_Data/L_565</SourceFilename>
       <SourceBand>1</SourceBand>
     </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="6">
     <SimpleSource>
       <SourceFilename relativeToVRT="1">HDF5:"COCTSOrigin.h5"://Geophysical_Data/L_670</SourceFilename>
       <SourceBand>1</SourceBand>
     </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="7">
     <SimpleSource>
       <SourceFilename relativeToVRT="1">HDF5:"COCTSOrigin.h5"://Geophysical_Data/L_750</SourceFilename>
       <SourceBand>1</SourceBand>
     </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="8">
     <SimpleSource>
       <SourceFilename relativeToVRT="1">HDF5:"COCTSOrigin.h5"://Geophysical_Data/L_865</SourceFilename>
       <SourceBand>1</SourceBand>
     </SimpleSource>
  </VRTRasterBand>
  
</VRTDataset>
```

因为是xml文件，本身就具有自描述性，这里也就不解释了。

构造完毕后，可以直接使用gdalwarp工具校正了：

.\gdalwarp.exe -geoloc -t\_srs \'+proj=aea +lat\_1=15 +lat\_2=65 +lat\_0=30 +lon\_0=95 +x\_0=0 +y\_0=0 +datum=WGS84 +units=m +no\_defs\' E:\testdata\COCTSOrigin.vrt E:\GEOREFCOCTS.tif