使用带有 Python 的 GDAL 裁剪光栅文件

Posted 2023-03-23

【中文标题】使用带有 Python 的 GDAL 裁剪光栅文件

【英文标题】：Cropping a raster file using GDAL w/ Python

【发布时间】：2018-07-09 21:48:46

【问题描述】：

编辑：部分解决方案

在命令行中使用 gdal_translate 似乎可以解决问题，即使 Python 绑定不起作用。

这有助于裁剪 GeoTiff，从顶部和左侧移除 300 像素的填充，并保留下一个 2000x2000 像素。 gdal_translate -srcwin 300 300 2000 2000 input.tif output.tif

原始问题

我花了很长时间试图弄清楚这一点。

我有一系列 GEOTiff 格式的卫星图像。每个图像都有一个 300 像素的缓冲区，它与旁边的图像重叠。

目标：我正在尝试从每个图像中裁剪出 300 像素的缓冲区，然后将其用作 GDAL 的光栅。

约束： 1）我需要保留与文件关联的所有元数据和坐标系统信息 2）我想完全用 Python 来做这个（没有命令行）

我尝试过的失败：

1)使用gdal.translate's srcWin函数：

raster_data = gdal.Open('image.tif')
x_size = raster_data.RasterXSize
y_size = raster_data.RasterYSize
raster_data_unpadded = gdal.Translate('temp.tif', raster_data,
                                       srcWin = [padding,padding, 
                                       x_size - padding*2, 
                                       y_size-padding*2])

问题：这会产生没有数据的黑色图像

2) 使用PIL 裁剪图像，然后另存为TIF

 from PIL import Image
 img = Image.open(image.tif)
 x_size, y_size = img.size
 box = (padding, padding, x_size-padding, y_size - padding)
 img_unpadded = img.crop(box)
 img_unpadded.save('unpadded_image.tif')

问题：PIL 无法保存文件。它只是挂起。尝试另存为“.tiff”文件会产生错误“编码器错误 -9”

3) 使用Rasterio

with rasterio.open("image.tif") as src:

    out_image, out_transform = mask(src, geoms, crop=True)

out_meta = src.meta.copy()

问题：Rasterio 似乎不接受像素格式的遮罩（例如 300 像素）。仅在文件的坐标系中采用几何图形，例如多边形。我不知道如何在像素和坐标之间进行转换。

【问题讨论】：

你能详细说明你到目前为止的尝试吗？如果它只有几行代码，您应该通过将其编辑到项目符号 3 下方来修改您的问题。

有什么理由蔑视命令行？与 python gdal 相反，它似乎正在工作

确认它使用 gdal 命令行工作。 gdal_translate -srcwin 300 300 2000 2000 input.tif output.tif [创建一个裁剪的 GeoTIFF 300 像素从顶部和左侧裁剪，高和长 2000 像素]。奇怪的是，同样的事情不适用于 Python 绑定......

【参考方案1】：

以下应该有效：

import rasterio
from rasterio.windows import Window

with rasterio.open('input.tif') as src:
    window = Window(padding, padding, src.width - 2 * padding, src.height - 2 * padding)

    kwargs = src.meta.copy()
    kwargs.update(
        'height': window.height,
        'width': window.width,
        'transform': rasterio.windows.transform(window, src.transform))

    with rasterio.open('cropped.tif', 'w', **kwargs) as dst:
        dst.write(src.read(window=window))

【参考方案2】：

这就是你通过gdal.Translate在python中使用gdal_translate的方式。

最佳选择：projwin

最简单的方法是使用projwin 标志，它有4 个值：

window = (upper_left_x, upper_left_y, lower_right_x, lower_right_y)

这些值在地图坐标中。输入文件的边界可以在命令行通过gdalinfo input_raster.tif获得。

注意：对于许多坐标系，ymax 实际上小于ymin，因此使用“upper_left”和“lower_right”来标识坐标很重要，而不是“max”和“分钟。”由于这种差异，Akin 的回答对我不起作用。

完整的解决方案是：

from osgeo import gdal

upper_left_x = 699934.584491
upper_left_y = 6169364.0093
lower_right_x = 700160.946739
lower_right_y = 6168703.00544
window = (upper_left_x,upper_left_y,lower_right_x,lower_right_y)

gdal.Translate('output_crop_raster.tif', 'input_raster.tif', projWin = window)

附加选项：srcwin

srcwin 是另一个 gdal_translate 标志，类似于projwin，但通过偏移量和大小接收像素和线窗口，而不是使用地图坐标边界。你会这样使用它，但 OP 似乎对这种方法有问题，所以最好只使用projwin：

window = (offset_x, offset_y, size_x, size_y)
gdal.Translate('output_crop_raster.tif', 'input_raster.tif', srcWin = window)