Suite

Raster d'écrêtage avec couche vectorielle à l'aide de GDAL

Raster d'écrêtage avec couche vectorielle à l'aide de GDAL


J'ai installé GDAL à l'aide du programme d'installation Osgeo. Comment puis-je découper une couche raster avec une couche vectorielle par programmation ? Existe-t-il une API GDAL disponible qui puisse m'aider? J'utilise Python.


Je ne suis pas sûr de l'API gdal, il y avoid* GDALWarpOptions::hCutlinedans les options Warp référencées dans le didacticiel de l'API Warp, mais pas d'exemples explicites. Êtes-vous sûr d'avoir besoin d'une réponse par programmation ? Les utilitaires de ligne de commande peuvent le faire directement :

  1. créer un fichier de formes contenant uniquement le polygone de découpage de la zone d'intérêt
  2. utilisation ogrinfo pour déterminer l'étendue du fichier de formes d'écrêtage
  3. utilisation gdal_translate pour découper à l'étendue de la forme
  4. utilisation gdalwarp avec-ligne de coupeparamètre

Les étapes 2 et 3 sont pour l'optimisation, vous pouvez vous en tirer avec justegdalwarp -ligne de coupe….

Voir Dé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 coupe peut être vu dans le didacticiel de Michael Corey sur la création d'ombrages pour Mapnik.


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

Le commentaire de Mike Toews est assez utile mais vous pourriez 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, c'est-à-dire aussi simple que possible et ne nécessitant pas de nombreuses dépendances externes. L'utilisation de Python Imaging Library 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 depuis (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 est la même), mais j'ai aimé l'idée de rendre le raster également plus petit (par exemple -crop_to_cutline), alors j'ai adoptémonde2Pixelde 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) # choisissez tous les points de données à l'intérieur des limites carrées de l'AOI, # remplacez tout autres points avec NULL craster.cdata= np.choose(np.flipud(craster.mask), (craster.data, -9999)) # redimensionner l'ensemble de données à la taille du polygone carré craster.ccdata=craster.cdata [ulY:lrY, ulX:lrX] craster.writer("ccdata2m.asc",craster.ccdata, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=Fa lse) # dans la deuxième étape, nous choisissons à nouveau tous les points de données qui sont à l'intérieur des # sommets englobants de l'AOI # devons redéfinir nos points raster craster.xllcorner, craster.yllcorner = minX, minY craster.xurcorner, craster.yurcorner = maxX , maxY craster.fillrasterpoints(10,10) craster.getmask(craster.boundingvertices) # juste un wrapper autour de 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) # écrire le raster sur le disque craster.writer("ccdata2m_clipped.asc", craster.data, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)

pour une description complète duclasse MaskRasteret ses méthodes, voir le github de mon projet.

En utilisant ce code, vous devrez toujours utiliser GDAL. Cependant, le plan est d'utiliser à l'avenir du Python pur là où je peux, car le public visé par 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…).


Joel Lawhead de GeospatialPython a un exemple complet de python dans Clip raster à l'aide de shapefile, un didacticiel bien écrit. Vous devrez installer la bibliothèque d'images Python (PIL) qui n'est pas incluse dans Osgeo4W (pour laquelle vous devrez peut-être ajouter o4w python au registre Windows pour que le programme d'installation fonctionne).


Voir la vidéo: 201 - Working with geotiff files using rasterio in python also quick demo of NDVI calculation