在进行遥感影像处理的时候，我们经常需要进行裁剪的工作，来看看如何使用GDAL工具进行这项操作吧！
参考资料：
GDAL: gdalwarp
GDAL: gdal_translate
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参数，如果进行矩形裁剪的话，只有GT(0)和GT(3)参数会有变化，即需要重新计算裁剪以后的左上角坐标即可。
下面给出使用Python对MODIS反射率的第一波段进行裁剪的代码：
from osgeo import gdal
import numpy as np
image_name = ('HDF4_EOS:EOS_GRID:'
'"MOD09GA.A2018349.h26v05.006.2018351030314.hdf":'
'MODIS_Grid_500m_2D:sur_refl_b01_1')
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
原文地址：
https://theonegis.github.io/geos/%E6%A0%85%E6%A0%BC%E6%95%B0%E6%8D%AE%E8%A3%81%E5%89%AA/index.html
在进行遥感影像处理的时候，我们经常需要进行裁剪的工作，来看看如何使用GDAL工具进行这项操作吧！参考资料：GDAL: gdalwarpGDAL: gdal_translateGDAL/OGR Python API使用GDAL命令GDAL提供了两个命令可以用于影像的裁剪：gdalwarp和gdal_translate，两个命令中我更推荐使用后者。gdalwarp命令可以使用-te制定裁剪范围。默认是在原数据的坐标系下的xmin ymin xmax ymax，当然我们也可以使用-te_srs参数指
本文是按照图像的像素坐标进行裁剪，需要自定义图像的初始像素点和需要裁剪的尺寸大小，并带有地理坐标信息。
#include <gdal.h>
#include <gdal_priv.h>
#include <iostream>
1、geotiff转换为GMT所需的grd文件
gdal_translate -of GMT -a_nodata NAN input.tif output.tif
2、有时候geotiff不止一个波段，需要转成单波段的geotiff才能进一步转化为GMT的grd文件
gdal_translate -b 1 input.tif output.tif
3、更改投影系统，比如把横轴墨卡托投影转换为带有经纬度的WGS84坐标系
gdalwrap -s_srs def -t_srs def input.tif out
本文收录了一些关于栅格的操作。GDAL对于影像的操作中，有很多需要注意的地方，可能由于影像的差异，所使用的的参数也存在一些差异，这里都进行了列举。希望能对大家有所帮助。
（1）获取影像信息
获取影像信息参数设置
public static string Gdal_Info(string filePath)
StringBuilder builder = new StringBuilder();
//builder.Append("-proj4 ");
//builder.App
利用shapefile对栅格图像进行裁剪.
; :Syntax
; RasterSubsetViaShapefile, Fid, shpFile=string, [pos=array],
; [inside={0|1}], [outFile=｛string|variable}], [r_fid=variable]
; :Params:
; Fid -- 输入文件FID
; 注：可通过ENVI_OPEN_FILE、ENVI_SELECT、ENVIRasterToFID等获取
; :Keywords:
; pos -- 保留波段索引数组（可选），默认保留所有波段。
; shpFile -- 用于裁剪的shapefile完整路径
; inside -- 保留shp文件外或内（可选，0或1），默认保留内部。
; 注：设置0时，保留外部；设置1时，保留内部。
; outFile -- 裁减结果文件路径（可选）
; 注：如果不设置或设置为变量，则裁剪结果保存在临时目录中，outFile将保存输出文件名
; 注：如果设置为文件路径，则裁剪结果保存在指定路径中
; r_fid -- 返回裁剪结果文件FID，如果范围-1，则表示裁剪失败。
本文讲解如何使用GDAL处理地理图像时，通过使用行列号计算和转换成tiff图像的地理坐标：
tif中坐标计算的方法如下，其中Col表示该坐标点处图像的列号，ROW表示该坐标点处图像的行号，
比如图像左上角Col为0，ROW为0，图像右下角Col为图像宽度，ROW为图像高度。
Xgeo = GT(0) + Col*GT(1) + Row*GT(2)
Ygeo = GT(3) + Col*GT(4) + Row*GT(5)
#include <gdal_pri..
使用gdal生成COG需要用到gdaladdo.exe与gdal_translate.exe参考网站是wiki的COG简介与部分操作：
https://trac.osgeo.org/gdal/wiki/CloudOptimizedGeoTIFF#NOTEgdaladdo.exe用来提前处理geotiff，生成概览层。这一步是生成COG的前提，因为COG内部得有概览层，而输入的geotiff不一定有。故在转换之前得有这一步。操作如下：
进入gdaladdo所在目录（C:\programs\gmt6\bin）。
# 定义输入矢量数据路径和栅格数据路径
vector_path = 'path/to/vector.shp'
raster_path = 'path/to/raster.tif'
# 打开矢量数据文件并获取几何信息
vector_ds = ogr.Open(vector_path)
layer = vector_ds.GetLayer()
feature = layer.GetFeature(0)
geometry = feature.GetGeometryRef()
# 打开栅格数据文件并获取地理参考和变换信息
raster_ds = gdal.Open(raster_path)
geo_transform = raster_ds.GetGeoTransform()
proj = osr.SpatialReference()
proj.ImportFromWkt(raster_ds.GetProjection())
# 将矢量数据的几何信息转换为栅格数据坐标系下的坐标
minX, maxX, minY, maxY = layer.GetExtent()
ulX, ulY = gdal.ApplyGeoTransform(geo_transform, minX, maxY)
lrX, lrY = gdal.ApplyGeoTransform(geo_transform, maxX, minY)
# 计算裁剪后的栅格数据的大小和地理参考
x_pixels = int((lrX - ulX) / geo_transform[1])
y_pixels = int((lrY - ulY) / geo_transform[5])
clip_proj = raster_ds.GetProjection()
# 创建输出栅格数据文件
driver = gdal.GetDriverByName('GTiff')
clip_raster_path = 'path/to/clip_raster.tif'
clip_raster_ds = driver.Create(clip_raster_path, x_pixels, y_pixels, 1, gdal.GDT_Float32)
clip_raster_ds.SetGeoTransform((ulX, geo_transform[1], 0, ulY, 0, geo_transform[5]))
clip_raster_ds.SetProjection(clip_proj)
# 裁剪栅格数据
gdal.Warp(clip_raster_ds, raster_ds, cutlineDSName=vector_path, cropToCutline=True)
# 关闭文件
clip_raster_ds = None
raster_ds = None
vector_ds = None
请注意，此代码假定输入矢量数据为多边形，并且只裁剪了栅格数据的第一个波段。如果需要裁剪多个波段，则需要使用适当的循环来处理每个波段，并将结果保存到多波段栅格数据中。此外，代码还需要更多的错误检查和边缘情况的处理。
GDAL中的地理坐标系、投影坐标系及其相互转换
fuuhoo:
GDAL中的地理坐标系、投影坐标系及其相互转换
陨星落云:
GDAL中的地理坐标系、投影坐标系及其相互转换
并行遥感:
