这个问题主要是关于使用Python、GDAL库、OGR库等工具，通过shapefile来对GTiff图像进行裁剪和掩膜。提问者想要利用shapefile中的多边形对大型GTiff图像进行裁剪，并将裁剪后的图像保存到指定目录中。同时，他希望裁剪后的图像具有正确的地理位置和值，并且不会出现所谓的“右阴影”问题。
2、解决方案
代码例子：
import os, sys, time, Tkinter as Tk, tkFileDialog
import operator
from osgeo import gdal, gdalnumeric, ogr, osr
import Image, ImageDraw
def SelectFile(req = 'Please select a file:', ft='txt'):
""" Customizable file-selection dialogue window, returns list() = [full path, root path, and filename]. """
try: # Try to select a csv dataset
foptions = dict(filetypes=[(ft+' file','*.'+ft)], defaultextension='.'+ft)
root = Tk.Tk(); root.withdraw(); fname = tkFileDialog.askopenfilename(title=req, **foptions); root.destroy()
return [fname]+list(os.path.split(fname))
except: print "Error: {0}".format(sys.exc_info()[1]); time.sleep(5); sys.exit()
def rnd(v, N): return int(round(v/float(N))*N)
def rnd2(v): return int(round(v))
# Raster image to clip
rname = SelectFile('Please select your TIF DEM:',ft='tif')
raster = rname[2]
print 'DEM:', raster
os.chdir(rname[1])
# Polygon shapefile used to clip
shp = SelectFile('Please select your shapefile of catchments (requires Name field):',ft='shp')[2]
print 'shp:', shp
qs = raw_input('Do you want to stretch the image? (y/n)')
qs = True if qs == 'y' else False
# Name of base clip raster file(s)
if not os.path.exists('./clip/'): os.mkdir('./clip/')
output = "/clip/clip"
#
def imageToArray(i):
"""
Converts a Python Imaging Library array to a
gdalnumeric image.
"""
a=gdalnumeric.fromstring(i.tostring(),'b')
a.shape=i.im.size[1], i.im.size[0]
return a
def arrayToImage(a):
"""
Converts a gdalnumeric array to a
Python Imaging Library Image.
"""
i=Image.fromstring('L',(a.shape[1],a.shape[0]), (a.astype('b')).tostring())
return i
def world2Pixel(geoMatrix, x, y, N= 5, r=True):
"""
Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
the pixel location of a geospatial coordinate
"""
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
rtnX = geoMatrix[2]
rtnY = geoMatrix[4]
if r:
pixel = int(round(x - ulX) / xDist)
line = int(round(ulY - y) / xDist)
else:
pixel = int(rnd(x - ulX, N) / xDist)
line = int(rnd(ulY - y, N) / xDist)
return (pixel, line)
def histogram(a, bins=range(0,256)):
"""
Histogram function for multi-dimensional array.
a = array
bins = range of numbers to match
"""
fa = a.flat
n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins)
n = gdalnumeric.concatenate([n, [len(fa)]])
hist = n[1:]-n[:-1]
return hist
def stretch(a):
"""
Performs a histogram stretch on a gdalnumeric array image.
"""
hist = histogram(a)
im = arrayToImage(a)
lut = []
for b in range(0, len(hist), 256):
# step size
step = reduce(operator.add, hist[b:b+256]) / 255
# create equalization lookup table
n = 0
for i in range(256):
lut.append(n / step)
n = n + hist[i+b]
im = im.point(lut)
return imageToArray(im)
#
srcImage = gdal.Open(raster)
geoTrans_src = srcImage.GetGeoTransform()
#print geoTrans_src
pxs = int(geoTrans_src[1])
srcband = srcImage.GetRasterBand(1)
ndv = -9999.0
#ndv = 0
# Create an OGR layer from a boundary shapefile
shapef = ogr.Open(shp)
lyr = shapef.GetLayer()
minXl, maxXl, minYl, maxYl = lyr.GetExtent()
ulXl, ulYl = world2Pixel(geoTrans_src, minXl, maxYl)
lrXl, lrYl = world2Pixel(geoTrans_src, maxXl, minYl)
#poly = lyr.GetNextFeature()
for poly in lyr:
pnm = poly.GetField("Name")
# Convert the layer extent to image pixel coordinates
geom = poly.GetGeometryRef()
#print geom.GetEnvelope()
minX, maxX, minY, maxY = geom.GetEnvelope()
geoTrans = geoTrans_src
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)
# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)
# Load the source data as a gdalnumeric array
#srcArray = gdalnumeric.LoadFile(raster)
clip = gdalnumeric.BandReadAsArray(srcband, xoff=ulX, yoff=ulY, win_xsize=pxWidth, win_ysize=pxHeight)
#clip = srcArray[:, ulY:lrY, ulX:lrX]
# Create a new geomatrix for the image
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY
# Map points to pixels for drawing the
# boundary on a blank 8-bit,
# black and white, mask image.
points = []
pixels = []
#geom = poly.GetGeometryRef()
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 = gdalnumeric.choose(mask, (clip, 0)).astype(gdalnumeric.uint8)
clip = gdalnumeric.choose(mask, (clip, ndv)).astype(gdalnumeric.numpy.float)
#
#for i in range(3):
# clip[i,:,:] = stretch(clip[i,:,:])
if qs: clip[:,:] = stretch(clip[:,:])
# Save ndvi as tiff
outputi = rname[1]+output+'_'+pnm+'.tif'
#gdalnumeric.SaveArray(clip, outputi, format="GTiff", prototype=srcImage)
driver = gdal.GetDriverByName('GTiff')
DataSet = driver.Create(outputi, pxWidth, pxHeight, 1, gdal.GDT_Float64)
#DataSet = driver.Create(outputi, pxWidth, pxHeight, 1, gdal.GDT_Int32)
DataSet.SetGeoTransform(geoTrans)
Projection = osr.SpatialReference()
