Raster de découpage avec couche vectorielle à l'aide de GDAL

26

J'ai installé GDAL à l'aide du programme d'installation d'Osgeo. Comment découper un calque raster avec un calque vectoriel par programmation? Existe-t-il une API GDAL disponible qui peut m'aider avec cela? J'utilise Python.

python gdal Vicky
la source

Réponses:

13

Je ne suis pas sûr de l'api gdal, il y a void* GDALWarpOptions::hCutlinedans les options Warp référencées dans le tutoriel Warp API , mais pas d'exemples explicites. Êtes-vous sûr d'avoir besoin d'une réponse programmatique? Les utilitaires de ligne de commande peuvent le faire hors de la boîte:

  1. créer un fichier de formes contenant uniquement la zone d'intérêt polygone de découpage
  2. utiliser ogrinfopour déterminer l'étendue du fichier de formes de découpage
  3. utiliser gdal_translatepour couper les étendues de forme
  4. utiliser gdalwarpavec -cutlineparamètre

Les étapes 2 et 3 sont pour l'optimisation, vous pouvez vous en tirer avec juste gdalwarp -cutline ....

Voir Couper des rasters avec GDAL en utilisant des polygones de Linfinity pour une solution basée sur Linux, le tout regroupé dans un seul script. Un autre exemple de ligne de démarcation peut être vu dans le tutoriel de Michael Corey créant des nuances pour Mapnik .

Matt Wilkie
la source
Matt, vous vous souvenez peut-être que trac.osgeo.org/gdal/ticket/1599 semble que la ligne de démarcation remplit cela
Mike T
10

Il semble que ce sujet revient toujours. Je ne savais pas moi-même que GDAL> 1.8 est si avancé qu'il vous offre déjà une gestion équitable de la ligne de commande pour effectuer cette tâche.

Le commentaire de Mike Toews est assez utile mais vous pouvez simplement faire par exemple:

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

Vous pouvez envelopper cette commande dans un script python avec l'excellent module de sous-processus .

Une chose qui était vraiment problématique pour moi est que je devais fournir une solution minimale à ce problème, ce qui signifie aussi simple que possible et ne nécessite pas de nombreuses dépendances externes. L'utilisation de la bibliothèque d'imagerie Python comme dans le tutoriel de Joel Lawhead est soignée, mais j'ai trouvé la solution suivante: utiliser des tableaux masqués Numpy.
Je ne sais pas si c'est mieux, mais c'est ce que je savais à l'époque (il y a 3 ans ...).
À l'origine, j'ai créé une zone de données valide à l'intérieur du raster d'origine (par exemple, l'étendue du raster en sortie où il était le même), mais j'ai aimé l'idée de rendre le raster également plus petit (par exemple, -crop_to_cutline), j'ai donc adopté world2Pixelde Joel Lawhead. Voici ma propre solution:

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)

pour une description complète des class MaskRasterméthodes et c'est, voir le github de mon projet .

En utilisant ce code, vous devrez toujours utiliser GDAL. Cependant, le plan est d'utiliser à l'avenir du pur Python où je peux, car l'audience visée de mon logiciel a des difficultés avec trop de dépendances (j'utilise Debian pour développer le logiciel, et les clients utilisent Windows 7 ...).

Oz123
la source
J'aime l'exemple de ligne de commande que vous avez donné, mais pouvez-vous expliquer ce que fait l'argument -crop_to_cutline? Je ne sais pas quel est son objectif étant donné que le fichier de formes de découpage est spécifié par -cutline.
hendra
1
l'option -cutline découpe le raster dans le cadre de délimitation interne de la couche de polygones. Par exemple, s'il est plus petit dans la mesure où le raster en sortie sera également plus petit. Sans cela, le raster en sortie sera de la même taille que l'original, mais avec NULL dans tous les points en dehors de la zone qui vous intéresse.
Oz123