版权声明:本文为博主原创文章,转载请注明原文出处!
写作时间:2019-03-22
在进行遥感影像处理的时候,我们经常需要进行裁剪的工作,来看看如何使用GDAL工具进行这项操作吧!
参考资料:
使用GDAL命令
GDAL提供了两个命令可以用于影像的裁剪:gdalwarp
和gdal_translate
,两个命令中我更推荐使用后者。
gdalwarp
命令可以使用-te
制定裁剪范围。默认是在原数据的坐标系下的xmin ymin xmax ymax
,当然我们也可以使用-te_srs
参数指定-te
参数所在的坐标系。
为什么不推荐gdalwarp
命令呢?这是因为gdalwarp
命令只提供了根据坐标系的范围进行裁剪,而不支持根据行列号的裁剪。这时候我们可以求助于gdal_translate
命令。
gdal_transalte
命令即支持使用-srcwin
参数指定行列号范围xoff yoff xsize ysize
,也支持使用-projwin
参数指定原数据坐标系下的范围ulx uly lrx lry
。同时提供参数-projwin_srs
可以用于指定-projwin
参数所在的坐标系,即跟gdalwarp
命令中的-te_srs
参数类似。
下面给出一个示例:
gdal_translate -of "GTiff" -srcwin 10 10 256 256 -a_scale 1 HDF4_EOS:EOS_GRID:"MOD09GA.A2018349.h26v05.006.2018351030314.hdf":MODIS_Grid_500m_2D:sur_refl_b01_1 sr_1.tif
这行命令将MODIS数据中的反射率的第一波段进行裁剪,起点为第10行第10列,输出大小为256$\times$255,输出格式为TIFF。
注意这行命令有一个-a_scale 1
参数,这个参数指定了裁剪过程不要对DN值进行缩放。如果不加这个值得话,输出图像的DN值会被根据原数据的Scale=10000
放大10000倍。
使用Python代码
对于使用Python代码进行裁剪,我们有两种方法:
- 第一就是对命令行对应的借口直接进行调用。这个最直接最简单。
- 第二就是首先自己选择出需要裁剪的区域,然后计算裁剪区域的GeoTransform的系数,最后将投影和GeoTransform系数赋值给裁剪子区域,写入输出文件。
我们知道GDAL中使用了六参数模型存储GeoTransform参数,如果进行矩形裁剪的话,只有GT(0)
和GT(3)
参数会有变化,即需要重新计算裁剪以后的左上角坐标即可。
下面给出使用Python对MODIS反射率的第一波段进行裁剪的代码:
1 | from osgeo import gdal |