OK, after messing around a bit, I changed the script from the site hyperlink in the second line of the comment. The purpose of the script is to pin / mask a LARGE raster (that is, which cannot fit into a 32-bit Python 2.7.5 application) in GTiff format with a shapefile with several polygons (each with an "Name" entry) and save the trimmed rasters in the “clip” subdirectory, where each masked grid is named after each polygon “Name”. Like the original script, it assumes that the GTiff and shapefile are in the same projection and overlap correctly and process ~ 100 masks / sec. However, I changed that the script to 1) works with a grid of heights with a floating point, 2) loads only the window of a large grid into memory, which is limited by the current polygon (i.e.Reduces the load on memory), 2) exports GTiff, which has the right (i.e. not shifted) geolocation and value.
HOWEVER, I am having problems with each net with a mask, which will be called a "right-handed shadow." That is, for each vertical line in the polygon, where the right side of this line is outside the specified polygon, the hidden grid will include one raster cell to the right of this polygonal side.
So my question is: what am I doing wrong, which gives the masked grid the right shadow?
I will try to figure out how to send the shapefile and tif example so that others can play. The code below also has comment lines for images with integer satellites (for example, as in the source code from geospatialpython.com).
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:
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))
rname = SelectFile('Please select your TIF DEM:',ft='tif')
raster = rname[2]
print 'DEM:', raster
os.chdir(rname[1])
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
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 = reduce(operator.add, hist[b:b+256]) / 255
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()
pxs = int(geoTrans_src[1])
srcband = srcImage.GetRasterBand(1)
ndv = -9999.0
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)
for poly in lyr:
pnm = poly.GetField("Name")
geom = poly.GetGeometryRef()
minX, maxX, minY, maxY = geom.GetEnvelope()
geoTrans = geoTrans_src
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)
clip = gdalnumeric.BandReadAsArray(srcband, xoff=ulX, yoff=ulY, win_xsize=pxWidth, win_ysize=pxHeight)
geoTrans = list(geoTrans)
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 = gdalnumeric.choose(mask, (clip, ndv)).astype(gdalnumeric.numpy.float)
if qs: clip[:,:] = stretch(clip[:,:])
outputi = rname[1]+output+'_'+pnm+'.tif'
driver = gdal.GetDriverByName('GTiff')
DataSet = driver.Create(outputi, pxWidth, pxHeight, 1, gdal.GDT_Float64)
DataSet.SetGeoTransform(geoTrans)
Projection = osr.SpatialReference()
Projection.ImportFromWkt(srcImage.GetProjectionRef())
DataSet.SetProjection(Projection.ExportToWkt())
DataSet.GetRasterBand(1).WriteArray(clip)
DataSet.GetRasterBand(1).SetNoDataValue(ndv)
print 'Saved:', outputi
del mask, clip, geom
del driver, DataSet
del shapef, srcImage, srcband