# 栅格数据裁剪

作者：阿振

邮箱：<tanzhenyugis@163.com>

博客：<https://blog.csdn.net/theonegis/article/details/80089375>

修改时间：2019-03-22

声明：本文为博主原创文章，转载请注明原文出处

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

参考资料：

1. [GDAL: gdalwarp](https://www.gdal.org/gdalwarp.html)
2. [GDAL: gdal\_translate](https://www.gdal.org/gdal_translate.html)
3. [GDAL/OGR Python API](https://gdal.org/python/)

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

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

```python
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
```