masque GTiff avec SHAPEFILE en python avec GDAL, OGR, etc

Question

OK, Après un peu de tripoter, je l'ai peaufiné un script de l'hyperlien du site dans la deuxième ligne de commentaire. Le but du script est de couper/masquer un raster LARGE (c'est-à-dire qui ne peut pas tenir dans une application Python 2.7.5 32 bits) au format GTiff avec un fichier de formes avec plusieurs polygones (chacun avec un enregistrement "Name") raster les rasters dans un sous-répertoire "clip", où chaque grille masquée porte le nom du "Nom" de chaque polygone. Comme le script original, il suppose que le GTiff et le shapefile sont dans la même projection et se chevauchent correctement, et il traite ~ 100 masks/sec. Cependant, j'ai modifié ce script pour 1) travailler avec une grille d'élévation à valeur flottante, 2) charger seulement la fenêtre de la plus grande grille dans la mémoire qui est délimitée par le polygone courant (ie réduire la charge de mémoire), 2) exporter GTiff's qui ont le droit (c'est-à-dire non décalé) géo-localisation et valeur.

Cependant, j'avoir des ennuis avec chaque grille masquée ayant un ce que je vais appeler une « ombre du côté droit ». C'est pour chaque ~ ligne verticale dans un polygone où le côté droit de cette ligne est en dehors du polygone donné, la grille masquée comprendra une cellule raster à la droite de ce polygone-côté.

Ainsi, ma question est, ce que je fais mal qui donne la grille masquée une ombre droite?

Je vais essayer de comprendre comment publier un exemple shapefile et TIF pour que les autres puissent se reproduire. Le code ci-dessous comporte également des lignes de commentaire pour des images satellites à valeurs entières (par exemple dans le code original de geospatialpython.com).

# RasterClipper.py - clip a geospatial image using a shapefile 
# http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html 
# http://gis.stackexchange.com/questions/57005/python-gdal-write-new-raster-using-projection-from-old 

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" 

# 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.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) 

# Also load as a gdal image to get geotransform 
# (world file) info 
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 the image using the mask 
    #clip = gdalnumeric.choose(mask, (clip, 0)).astype(gdalnumeric.uint8) 
    clip = gdalnumeric.choose(mask, (clip, ndv)).astype(gdalnumeric.numpy.float) 

    # This image has 3 bands so we stretch each one to make them 
    # visually brighter 
    #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() 
    Projection.ImportFromWkt(srcImage.GetProjectionRef()) 
    DataSet.SetProjection(Projection.ExportToWkt()) 
    # Write the array 
    DataSet.GetRasterBand(1).WriteArray(clip) 
    DataSet.GetRasterBand(1).SetNoDataValue(ndv) 

    # Save ndvi as an 8-bit jpeg for an easy, quick preview 
    #clip = clip.astype(gdalnumeric.uint8) 
    #gdalnumeric.SaveArray(clip, rname[1]+outputi+'.jpg', format="JPEG") 
    #print '\t\tSaved:', outputi, '-.tif, -.jpg' 
    print 'Saved:', outputi 
    del mask, clip, geom 
    del driver, DataSet 

del shapef, srcImage, srcband 

Answer 1

Cette fonctionnalité est déjà incorporée dans les utilitaires de ligne de commande gdal. Compte tenu de votre cas, je ne vois aucune raison pour laquelle vous voulez le faire vous-même en Python.

Vous pouvez passer en boucle les geomerties avec OGR et pour chacun appel gdalwarp avec les paramètres appropriés.

import ogr 
import subprocess 

inraster = 'NE1_HR_LC_SR_W_DR\NE1_HR_LC_SR_W_DR.tif' 
inshape = '110m_cultural\ne_110m_admin_0_countries_lakes.shp' 

ds = ogr.Open(inshape) 
lyr = ds.GetLayer(0) 

lyr.ResetReading() 
ft = lyr.GetNextFeature() 

while ft: 

    country_name = ft.GetFieldAsString('admin') 

    outraster = inraster.replace('.tif', '_%s.tif' % country_name.replace(' ', '_'))  
    subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, 
        '-crop_to_cutline', '-cwhere', "'admin'='%s'" % country_name]) 

    ft = lyr.GetNextFeature() 

ds = None 

J'ai utilisé quelques exemples de données de la Terre naturelle dans l'exemple ci-dessus, pour le Brésil la découpe ressemble:

Si vous ne souhaitez que pour recadrer l'image dans la zone du polygone et ne masque rien à l'extérieur, vous pouvez transformer le Shapefile afin qu'il contienne l'enveloppe des polygones. Ou simplement desserrer le fichier de formes et appelez gdal_translate avec -projwin pour spécifier la zone d'intérêt.