#1 楼
我不确定gdal api，在Warp API教程中引用的Warp Options中有void* GDALWarpOptions::hCutline，但没有明确的示例。您确定需要程序化答案吗？命令行实用程序可以立即使用：创建一个仅包含感兴趣区域裁剪多边形的shapefile
使用
ogrinfo确定裁剪shapefile的范围使用
gdal_translate来裁剪形状范围将
gdalwarp与-cutline参数一起使用步骤2和3是为了进行优化，仅需
gdalwarp -cutline ...即可完成。 有关基于Linux的解决方案，请参阅使用Linfinity的多边形使用GDAL剪切栅格，这些栅格都封装在一个脚本中。迈克尔·科里（Michael Corey）的教程为Mapnik创建山体阴影时，可以看到另一个插队示例。
#2 楼
似乎这个话题总是回来的。我本人并不知道GDAL> 1.8是否如此先进，它已经为您提供了公平的命令行处理能力来完成该任务。Mike Toews的评论非常有用，但是您可以例如执行以下操作：
gdalwarp -of GTiff -cutline DATA/area_of_interest.shp -cl area_of_interest -crop_to_cutline DATA/PCE_in_gw.asc data_masked7.tiff
您可以将此命令包装在python脚本中具有出色的子流程模块。
对我来说真正有问题的一件事是，我需要为该问题提供最小的解决方案，这意味着尽可能简单，并且不需要很多外部依赖。 Joel Lawhead的教程中的Python Imaging Library的使用非常简洁，但是我想出了以下解决方案：使用Numpy掩码数组。
我不知道它是否更好，但这是我所知道的（3年前...）。
我最初在原始栅格内创建了一个有效的数据区域（例如，输出栅格的范围相同），但是我喜欢将栅格也较小（例如-crop_to_cutline），因此我采用了Joel Lawhead的
world2Pixel。 这是我自己的解决方案：
def RasterClipper():
craster = MaskRaster()
contraster2 = 'PCE_in_gw.aux'
craster.reader("DATA/"+contraster2.replace('aux','asc'))
xres, yres = craster.extent[1], craster.extent[1]
craster.fillrasterpoints(xres, yres)
craster.getareaofinterest("DATA/area_of_interest.shp")
minX, maxX=craster.new_extent [0]-5,craster.new_extent[1]+5
minY, maxY= craster.new_extent [2]-5,craster.new_extent[3]+5
ulX, ulY=world2Pixel(craster.extent, minX, maxY)
lrX, lrY=world2Pixel(craster.extent, maxX, minY)
craster.getmask(craster.corners)
craster.mask=np.logical_not(craster.mask)
craster.mask.resize(craster.Yrange.size,craster.Xrange.size)
# choose all data points inside the square boundaries of the AOI,
# replace all other points with NULL
craster.cdata= np.choose(np.flipud(craster.mask), (craster.data, -9999))
# resise the data set to be the size of the squared polygon
craster.ccdata=craster.cdata[ulY:lrY, ulX:lrX]
craster.writer("ccdata2m.asc",craster.ccdata, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)
# in second step we rechoose all the data points which are inside the
# bounding vertices of AOI
# need to re-define our raster points
craster.xllcorner, craster.yllcorner = minX, minY
craster.xurcorner, craster.yurcorner = maxX, maxY
craster.fillrasterpoints(10,10)
craster.getmask(craster.boundingvertices) # just a wrapper around matplotlib.nxutils.points_in_poly
craster.data=craster.ccdata
craster.clip2(new_extent_polygon=craster.boundingvertices)
craster.data = np.ma.MaskedArray(craster.data, mask=craster.mask)
craster.data = np.ma.filled(craster.data, fill_value=-9999)
# write the raster to disk
craster.writer("ccdata2m_clipped.asc",craster.data, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)
有关
class MaskRaster及其方法的完整说明，请参见我项目的github。使用此代码，您仍然需要使用GDAL。但是，计划是将来在可能的地方使用纯Python，因为我的软件的目标受众在依赖项过多时会遇到困难（我使用Debian开发软件，而客户端使用Windows 7 ...）。
评论
我喜欢您提供的命令行示例，但是可以解释一下-crop_to_cutline参数的作用吗？我不确定剪切形状文件是由-cutline指定的。
– hendra
2012年7月23日在2:18
-cutline选项将栅格裁剪到多边形层的内部边界框。例如，如果范围较小，则输出栅格也将较小。如果没有此设置，则输出栅格的大小将与原始栅格大小相同，但是在您感兴趣的区域之外的所有点都将为NULL。
–Oz123
2012年7月23日在9:56
评论
马特，您可能还记得trac.osgeo.org/gdal/ticket/1599
– Mike T
2012年1月26日在1:59