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

GDAL使用(五)

admin3个月前 (02-14)代码相关158

在之前说的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:\testdata\MOD021KM.A2016352.0250.061.2017329005547.hdf

modisinfo.jpg

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

选取EV_1KM_RefSB数据集进行校正。在powershell中输入以下内容:

.\gdalwarp.exe -geoloc  -t_srs '+proj=longlat +datum=WGS84 +no_defs 'HDF4_EOS:EOS_SWATH:"E:\testdata\MOD021KM.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

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

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

COCTSinfo.jpg

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

数据路径:以红光波段为例:HDF5:"E:\testdata\COCTSOrigin.h5"://Geophysical_Data/L_670

纬度路径:HDF5:"E:\testdata\COCTSOrigin.h5"://Navigation_Data/Latitude

经度路径:HDF5:"E:\testdata\COCTSOrigin.h5"://Navigation_Data/Longitude

然后构造vrt文件:

<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

版权声明:本文由cyhu's essay发布,如需转载请注明出处。

本文链接:https://lovedm.club/?id=120

相关文章

C、C++、C#交换变量

C、C++、C#交换变量

最近重新看了看C和C++,觉得有些地方挺有意思。作为一开始不管什么资料都会用来做例子的一个程序,交换变量。不管在哪,常用的int,float,double类型的变量都是值类型的,作为参数传到函数(方法...

C# 委托

C# 委托

Action和Func是.NET类库的内置委托,以简便使用。Func有17个重载还可以使用delegate关键字创建委托下面的代码展示了这三种使用委托的方法namespace _202003...

C# Lambda表达式

简单用法,一句一句来,便于理解 Func<int, int, int> func = new Func<int, int, int>((int a, int b) =...

数值微分法(DDA)绘制直线

数值微分法(DDA)绘制直线

数值微分法(Digital Differential Analyzer)直接从直线的微分方程生成直线。详细的原理见以下链接:https://blog.csdn.net/weixin_43751983/...

偶然想到的一个问题。。。

偶然想到的一个问题。。。

今天突然想C#中,用数组中的Max()方法和直接通过比较获取最大值的时间谁快,于是试了试:       static v...

C语言 rename

在<stdio.h>头文件下以下内容来自:http://www.cplusplus.com/reference/cstdio/rename/?kw=renamerenameint ...