Ausschnittsraster mit Vektorebene unter Verwendung von GDAL

26

Ich habe GDAL mit dem Osgeo-Installationsprogramm installiert. Wie kann ich programmgesteuert eine Rasterebene mit einer Vektorebene beschneiden? Gibt es eine GDAL-API, die mir dabei helfen kann? Ich benutze Python.

Vicky
quelle

Antworten:

13

Ich bin mir nicht sicher, was die GDAL-API angeht, da void* GDALWarpOptions::hCutlinein den Warp-Optionen auf das Warp-API-Tutorial verwiesen wird , aber keine expliziten Beispiele. Sind Sie sicher, dass Sie eine programmatische Antwort benötigen? Die Befehlszeilen-Dienstprogramme können dies sofort tun:

  1. Erstellen Sie ein Shapefile, das nur das Beschneidungspolygon des gewünschten Bereichs enthält
  2. Verwenden Sie ogrinfodiese Option, um das Ausmaß des Clipping-Shapefiles zu bestimmen
  3. Verwenden Sie gdal_translatediese Taste, um die Ausmaße der Form anzupassen
  4. Verwenden Sie gdalwarpmit -cutlineParameter

Die Schritte 2 und 3 dienen der Optimierung gdalwarp -cutline ....

Weitere Informationen finden Sie unter Ausschneiden von Rastern mit GDAL mithilfe von Polygonen von Linfinity für linuxbasierte Lösungen, die alle in einem Skript zusammengefasst sind. Ein weiteres Beispiel für eine Schnittlinie ist in Michael Coreys Tutorial zu sehen, in dem Sie für Mapnik Schattierungen erstellen .

Matt Wilkie
quelle
Matt, Sie können sich vielleicht erinnern, dass trac.osgeo.org/gdal/ticket/1599 so aussieht, als ob die Cutline dies erfüllt
Mike T
10

Es scheint, dass dieses Thema immer wieder auftaucht. Ich selbst wusste nicht, dass GDAL> 1.8 so weit fortgeschritten ist, dass Sie für diese Aufgabe bereits eine faire Befehlszeilensteuerung haben.

Der Kommentar von Mike Toews ist ziemlich nützlich, aber Sie könnten einfach zum Beispiel tun:

gdalwarp -of GTiff -cutline DATA/area_of_interest.shp -cl area_of_interest  -crop_to_cutline DATA/PCE_in_gw.asc  data_masked7.tiff 

Sie können diesen Befehl in ein Python-Skript mit dem hervorragenden Unterprozessmodul einbinden.

Eine Sache, die für mich wirklich problematisch war, war, dass ich eine minimale Lösung für dieses Problem liefern musste, was bedeutet, dass es so einfach wie möglich ist und nicht viele externe Abhängigkeiten erfordert. Die Verwendung der Python-Imaging-Bibliothek wie im Tutorial von Joel Lawhead ist ordentlich, aber ich habe die folgende Lösung gefunden: Verwenden von Numpy-maskierten Arrays.
Ich weiß nicht, ob es besser ist, aber das war es, was ich damals wusste (vor 3 Jahren ...).
Ursprünglich habe ich einen gültigen Datenbereich innerhalb des ursprünglichen Rasters erstellt (z. B. das Ausmaß des Ausgabe-Rasters, in dem es dasselbe ist), aber ich mochte die Idee, das Raster auch kleiner zu machen (z. B. -crop_to_cutline), also habe ich world2Pixelvon Joel Lawhead übernommen. Hier ist meine eigene Lösung:

def RasterClipper():
    craster = MaskRaster()
    contraster2 = 'PCE_in_gw.aux'
    craster.reader("DATA/"+contraster2.replace('aux','asc'))
    xres, yres = craster.extent[1], craster.extent[1]
    craster.fillrasterpoints(xres, yres)
    craster.getareaofinterest("DATA/area_of_interest.shp")
    minX, maxX=craster.new_extent [0]-5,craster.new_extent[1]+5
    minY, maxY= craster.new_extent [2]-5,craster.new_extent[3]+5
    ulX, ulY=world2Pixel(craster.extent, minX, maxY)
    lrX, lrY=world2Pixel(craster.extent, maxX, minY)
    craster.getmask(craster.corners)
    craster.mask=np.logical_not(craster.mask)
    craster.mask.resize(craster.Yrange.size,craster.Xrange.size)
    # choose all data points inside the square boundaries of the AOI,
    # replace all other points with NULL
    craster.cdata= np.choose(np.flipud(craster.mask), (craster.data, -9999))
    # resise the data set to be the size of the squared polygon
    craster.ccdata=craster.cdata[ulY:lrY, ulX:lrX]
    craster.writer("ccdata2m.asc",craster.ccdata, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)
    # in second step we rechoose all the data points which are inside the
    # bounding vertices of AOI
    # need to re-define our raster points
    craster.xllcorner, craster.yllcorner = minX, minY
    craster.xurcorner, craster.yurcorner = maxX, maxY
    craster.fillrasterpoints(10,10)
    craster.getmask(craster.boundingvertices) # just a wrapper around matplotlib.nxutils.points_in_poly
    craster.data=craster.ccdata
    craster.clip2(new_extent_polygon=craster.boundingvertices)
    craster.data = np.ma.MaskedArray(craster.data, mask=craster.mask)
    craster.data = np.ma.filled(craster.data, fill_value=-9999)
    # write the raster to disk
    craster.writer("ccdata2m_clipped.asc",craster.data, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)

Eine vollständige Beschreibung der class MaskRasterund ihrer Methoden finden Sie im Github meines Projekts .

Mit diesem Code müssen Sie weiterhin GDAL verwenden. Der Plan ist jedoch, in Zukunft reines Python zu verwenden, wo ich kann, da die Zielgruppe meiner Software Schwierigkeiten mit zu vielen Abhängigkeiten hat (ich benutze Debian, um die Software zu entwickeln, und die Clients verwenden Windows 7 ...).

Oz123
quelle
Ich mag das Befehlszeilenbeispiel, das Sie gegeben haben, aber können Sie erklären, was das Argument -crop_to_cutline bewirkt? Ich bin mir nicht sicher, wozu es dient. Das Clipping-Shapefile wird durch -cutline angegeben.
Hendra
1
Mit der Option -cutline wird das Raster auf den inneren Begrenzungsrahmen der Polygonebene gekürzt. Wenn es beispielsweise in den Bereichen kleiner ist, ist auch das Ausgabe-Raster kleiner. Ohne diese Option hat das Ausgabe-Raster die gleiche Größe wie das Original, jedoch an allen Stellen außerhalb des Bereichs Ihres Interesses den Wert NULL.
Oz123