March 22, 2019

栅格数据裁剪

版权声明:本文为博主原创文章,转载请注明原文出处!
写作时间:2019-03-22

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

参考资料:

  1. GDAL: gdalwarp
  2. GDAL: gdal_translate
  3. GDAL/OGR Python API

使用GDAL命令

GDAL提供了两个命令可以用于影像的裁剪:gdalwarpgdal_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代码进行裁剪,我们有两种方法:

我们知道GDAL中使用了六参数模型存储GeoTransform参数,如果进行矩形裁剪的话,只有GT(0)GT(3)参数会有变化,即需要重新计算裁剪以后的左上角坐标即可。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
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