使用gdal python 接口根据多边形裁切数据

仅作记录, 因为网上找得没有一个能用得, 开发环境vs,  python 3.6, gdal3.1.4

from osgeo import gdal, gdalnumeric, ogr
from PIL import Image, ImageDraw
from osgeo import gdal_array
import os
#保存tif文件函数
import gdal
import numpy as np
gdal.UseExceptions()
 

# This function will convert the rasterized clipper shapefile
# to a mask for use within GDAL.
def imageToArray(i):
  """
  Converts a Python Imaging Library array to a
  gdalnumeric image.
  """
  a=gdalnumeric.fromstring(i.tobytes(),'b')
  a.shape=i.im.size[1], i.im.size[0]
  return a

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)
 
#
# EDIT: this is basically an overloaded
# version of the gdal_array.OpenArray passing in xoff, yoff explicitly
# so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
  ds =gdal_array.OpenArray(array)
 
  if ds is not None and prototype_ds is not None:
    if type(prototype_ds).__name__ == 'str':
      prototype_ds = gdal.Open( prototype_ds )
    if prototype_ds is not None:
      gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
  return ds


 
def main( shapefile_path, raster_path,outpath_path ):

  gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
  #gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
  # Also load as a gdal image to get geotransform
  # (world file) info
  srcImage = gdal.Open(raster_path)
  srcgeoTrans = srcImage.GetGeoTransform()
 
  # Create an OGR layer from a boundary shapefile
  shapef = ogr.Open(shapefile_path,1)
  lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
  lyr.ResetReading()
  poly = lyr.GetNextFeature()

  while poly:
        
      geom = poly.GetGeometryRef()
      minX, maxX, minY, maxY = geom.GetEnvelope()
      ulX, ulY = world2Pixel(srcgeoTrans, minX, maxY)
      lrX, lrY = world2Pixel(srcgeoTrans, maxX, minY)
      if ulX < 0:
          ulX =0
      if ulY < 0:
          ulY = 0
      if lrX < 0:
          lrX =0
      if lrY < 0:
          lrY = 0
      pxWidth = int(lrX - ulX)
      pxHeight = int(lrY - ulY)

      clip = gdalnumeric.LoadFile(raster_path,ulX,ulY,pxWidth,pxHeight)

      xoffset = ulX
      yoffset = ulY
      print ("Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset ))
 
      # Create a new geomatrix for the image
      geoTrans = list(srcgeoTrans)
      geoTrans[0] = minX
      geoTrans[3] = maxY
 

      points = []
      pixels = []
      pts = geom.GetGeometryRef(0)
      for p in range(pts.GetPointCount()):
       points.append((pts.GetX(p), pts.GetY(p)))
      for p in points:
       pixels.append(world2Pixel(geoTrans, p[0], p[1]))
      rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
      rasterize = ImageDraw.Draw(rasterPoly)
      rasterize.polygon(pixels, 0)
      mask = imageToArray(rasterPoly)
 
      # Clip the image using the mask
      clip = gdalnumeric.choose(mask, \
        (clip, 0)).astype(gdalnumeric.uint8)

 
      clip = clip.astype(gdalnumeric.uint8)
      out = outpath_path +"\\"+ str(poly.GetFID())+".jpg"
      gdalnumeric.SaveArray(clip, out, format="JPEG")
      field_index = lyr.GetLayerDefn().GetFieldIndex("Code")
      poly.SetField(field_index, "111")
      lyr.SetFeature(poly)
      poly = lyr.GetNextFeature()
  shapef.SyncToDisk()
  shapef.Destroy()
  gdal.ErrorReset()

 
if __name__ == '__main__': 
  #shapefile_path, raster_path 
  shapefile_path = r'G:\Polygon.shp' 
  raster_path = r'G:\a.img' 
  outpath_path = r'G:\21'
  main( shapefile_path, raster_path, outpath_path )