python图像切割成多边形_Python：用于剪切带有多边形的光栅图像的函数的代码检查...

我正在请求审核以下代码.我有一个空间参考图像和一个多边形.我编写了一个代码(见下文)来剪辑此图像以保存新图像(剪裁区域).此功能基于要素类的几何图形剪切栅格.基于几何体的剪切意味着您将使用要素类中所有要素的边界来剪切栅格,而不是这些要素的最小边界矩形

输入：多边形图层和一个或多个栅格图层

输出：新的栅格图层,剪切为多边形边界

import osgeo.gdal

import shapefile

import struct,numpy,pylab

import numpy as np

import ogr

import osr,gdal

from shapely.geometry import Polygon

import osgeo.gdal as gdal

import sys

from osgeo import gdal,gdalnumeric,ogr,osr

import Image,ImageDraw

def world2Pixel(geoMatrix,x,y):

"""

Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate

the pixel location of a geospatial coordinate

(http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html)

geoMatrix

[0] = top left x (x Origin)

[1] = w-e pixel resolution (pixel Width)

[2] = rotation,0 if image is "north up"

[3] = top left y (y Origin)

[4] = rotation,0 if image is "north up"

[5] = n-s pixel resolution (pixel Height)

"""

ulX = geoMatrix[0]

ulY = geoMatrix[3]

xDist = geoMatrix[1]

yDist = geoMatrix[5]

rtnX = geoMatrix[2]

rtnY = geoMatrix[4]

pixel = np.round((x - ulX) / xDist).astype(np.int)

line = np.round((ulY - y) / xDist).astype(np.int)

return (pixel,line)

def Pixel2world(geoMatrix,y):

ulX = geoMatrix[0]

ulY = geoMatrix[3]

xDist = geoMatrix[1]

yDist = geoMatrix[5]

coorX = (ulX + (x * xDist))

coorY = (ulY + (y * yDist))

return (coorX,coorY)

def RASTERClipByPolygon(inFile,poly,outFile):

# Open the image as a read only image

ds = osgeo.gdal.Open(inFile,gdal.GA_ReadOnly)

# Check the ds (=dataset) has been successfully open

# otherwise exit the script with an error message.

if ds is None:

raise SystemExit("The raster could not openned")

# Get image georeferencing information.

geoMatrix = ds.GetGeoTransform()

ulX = geoMatrix[0]

ulY = geoMatrix[3]

xDist = geoMatrix[1]

yDist = geoMatrix[5]

rtnX = geoMatrix[2]

rtnY = geoMatrix[4]

# get the WKT (= Well-known text)

dsWKT = ds.GetProjectionRef()

# get driver information

DriverName = ds.GetDriver().ShortName

# open shapefile (= border of are of interest)

shp = osgeo.ogr.Open(poly)

if len(shp.GetLayer()) != 1:

raise SystemExit('The shapefile must have exactly one layer')

# Create an OGR layer from a boundary shapefile

layer = shp.GetLayer(0)

feature = layer.GetNextFeature()

geometry = feature.GetGeometryRef()

# Make sure that it is a polygon

if geometry.GetGeometryType() != osgeo.ogr.wkbPolygon:

raise SystemExit('This module can only load polygon')

# get Extent of the clip area

X_min,X_max,Y_min,Y_max = layer.GetExtent()

# Convert the layer extent to image pixel coordinates

uldX,uldY = world2Pixel(geoMatrix,X_min,Y_max)

lrdX,lrdY = world2Pixel(geoMatrix,Y_min)

# Calculate the pixel size of the new image

pxWidth = int(lrdX - uldX)

pxHeight = int(lrdY - uldY)

# get the Coodinate of left-up vertex of the pixel

X_minPixel,Y_maxPixel = Pixel2world(geoMatrix,uldX,uldY)

# get polygon's vertices

pts = geometry.GetGeometryRef(0)

points = []

for p in range(pts.GetPointCount()):

points.append((pts.GetX(p),pts.GetY(p)))

pnts = np.array(points).transpose()

# work band by band

nBands = ds.RasterCount

# panchromatic

if nBands == 1:

band = ds.GetRasterBand(1)

# get nodata value

nodata = band.GetNoDataValue()

# convert band in Array

bandArray = band.ReadAsArray()

del band

# clip arrey

bandArray_Area = bandArray[uldY:lrdY,uldX:lrdX]

del bandArray

# Create 2D Polygon Mask. Mode 'L',not '1',because

# Numpy-1.5.0 / PIL-1.1.7 does not support the numpy.array(img)

# conversion nicely for bivalue images.

img = Image.new('L',(pxWidth,pxHeight),0)

target_ds = gdal.GetDriverByName(DriverName).Create(outFile,pxWidth,pxHeight,nBands,ds.GetRasterBand(1).DataType)

target_ds.SetGeoTransform((X_minPixel,xDist,rtnX,Y_maxPixel,rtnY,yDist))

pixels,line = world2Pixel(target_ds.GetGeoTransform(),pnts[0],pnts[1])

listdata = [(pixels[i],line[i]) for i in xrange(len(pixels))]

ImageDraw.Draw(img).polygon(listdata,outline=1,fill=1)

mask = numpy.array(img)

bandArray_Masked = bandArray_Area*mask

del bandArray_Area,mask

target_ds.GetRasterBand(nBands).WriteArray(bandArray_Masked)

target_ds.GetRasterBand(nBands).SetNoDataValue(nodata)

else:

img = Image.new('L',ds.GetRasterBand(1).DataType)

target_ds.SetGeoTransform((X_min,Y_max,fill=1)

mask = numpy.array(img)

for bandno in range(1,nBands+1):

band = ds.GetRasterBand(bandno)

nodata = band.GetNoDataValue()

# convert band in Array

bandArray = band.ReadAsArray()

del band

# clip arrey

bandArray_Area = bandArray[ulY:lrY,ulX:lrX]

del bandArray

bandArray_Masked = bandArray_Area*mask

target_ds.GetRasterBand(bandno).WriteArray(bandArray_Masked)

del bandArray_Area

target_ds.GetRasterBand(bandno).SetNoDataValue(nodata)

# set the reference info

if len(dsWKT) is 0:

# Source has no projection (needs GDAL >= 1.7.0 to work)

target_ds.SetProjection('LOCAL_CS["arbitrary"]')

else:

# Make the target raster have the same projection as the source

target_ds.SetProjection(dsWKT)

target_ds = None