栅格数据裁剪

缺乏、安全感 2022-03-01 16:14 342阅读 0赞

### 文章目录 ###

*   *  使用GDAL命令
     *  使用Python代码

写作时间：2019-03-22

在进行遥感影像处理的时候，我们经常需要进行裁剪的工作，来看看如何使用GDAL工具进行这项操作吧！

参考资料：

1.  [GDAL: gdalwarp][GDAL_ gdalwarp]
2.  [GDAL: gdal\_translate][GDAL_ gdal_translate]
3.  [GDAL/OGR Python API][GDAL_OGR Python API]

## 使用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参数][GDAL_GeoTransform]，如果进行矩形裁剪的话，只有`GT(0)`和`GT(3)`参数会有变化，即需要重新计算裁剪以后的左上角坐标即可。

下面给出使用Python对MODIS反射率的第一波段进行裁剪的代码：

from osgeo import gdal
    import numpy as np
    
    # API参考：https://gdal.org/python/
    # GDAL命令行参考：https://www.gdal.org/gdal_translate.html
    image_name = ('HDF4_EOS:EOS_GRID:'
                  '"MOD09GA.A2018349.h26v05.006.2018351030314.hdf":'
                  'MODIS_Grid_500m_2D:sur_refl_b01_1')
    
    # 第一种方式，也是最简单的方法：直接使用GDAL命令行对应的Python方法
    src: gdal.Dataset = gdal.Open(image_name)
    src = gdal.Translate('cropped_with_translate.tif', src, srcWin=[10, 10, 256, 256],
                         options=['-a_scale', '1'])
    del src
    
    # 第二种方式，自己选择出需要的像素，然后自己确定裁剪以后的空间参考关系，并写入到输出文件
    src: gdal.Dataset = gdal.Open(image_name)
    band: gdal.Band = src.GetRasterBand(1)
    subset: np.ndarray = band.ReadAsArray(10, 10, 256, 256)
    
    driver: gdal.Driver = gdal.GetDriverByName('GTiff')
    dst: gdal.Dataset = driver.Create('cropped_from_scratch.tif', 256, 256, 1, gdal.GDT_Int16)
    dst.SetProjection(src.GetProjection())
    trans = list(src.GetGeoTransform())
    trans[0] -= -10 * trans[1]
    trans[3] -= -10 * trans[5]
    dst.SetGeoTransform(tuple(trans))
    
    band: gdal.Band = dst.GetRasterBand(1)
    band.WriteArray(subset)
    band.FlushCache()
    del src
    del dst

[GDAL_ gdalwarp]: https://www.gdal.org/gdalwarp.html
[GDAL_ gdal_translate]: https://www.gdal.org/gdal_translate.html
[GDAL_OGR Python API]: https://gdal.org/python/
[GDAL_GeoTransform]: https://blog.csdn.net/theonegis/article/details/80304873