Recorte de ráster con capa vectorial utilizando GDAL

24

He instalado GDAL utilizando el instalador de Osgeo. ¿Cómo puedo recortar una capa ráster con una capa vectorial mediante programación? ¿Hay alguna API de GDAL disponible que pueda ayudarme con esto? Estoy usando Python.

    
pregunta Vicky 08.11.2011 - 14:12

3 respuestas

11

No estoy seguro de la gi api, hay void* GDALWarpOptions::hCutline en las Opciones de deformación a las que se hace referencia desde < a href="http://www.gdal.org/warptut.html"> Tutorial de Warp API , pero no hay ejemplos explícitos. ¿Estás seguro de que necesitas una respuesta programática? Las utilidades de línea de comando pueden hacerlo fuera de la caja:

  1. cree un shapefile que contenga solo el área del polígono de recorte de interés
  2. use ogrinfo para determinar la extensión del archivo de forma de recorte
  3. use gdal_translate para enganchar a las extensiones de forma
  4. use gdalwarp con el parámetro -cutline

Pasos 2 y amp; 3 son para optimización, usted podría salir adelante con solo gdalwarp -cutline ... .

Vea Recorte de rasters con GDAL usando polígonos de Linfinity para Linux Solución todo envuelto en un script. Otro ejemplo de línea de corte se puede ver en el tutorial de Michael Corey creando sombreados para Mapnik .

    
respondido por el matt wilkie 08.11.2011 - 19:42
9

Joel Lawhead de GeospatialPython tiene un ejemplo completo de python en Trama de clips usando shapefile , un tutorial bien escrito. Deberá instalar la Biblioteca de imágenes de Python (PIL) que no está incluida en Osgeo4W (para que podría necesitar agregar o4w python al registro de Windows para que el programa de instalación funcione).

    
respondido por el matt wilkie 08.11.2011 - 19:50
9

Parece que este tema siempre está volviendo. Yo mismo no sabía que GDAL > 1.8 es tan avanzado que ya le da un manejo justo en la línea de comandos para realizar esa tarea.

El comentario de Mike Toews es bastante útil, pero puedes hacerlo, por ejemplo:

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

Puede envolver este comando dentro de un script de Python con el excelente subprocess .

Una cosa que fue realmente problemática para mí es que necesitaba proporcionar una solución mínima a ese problema, es decir, lo más simple posible y no requiere muchas dependencias externas. El uso de Python Imaging Library como en El tutorial de Joel Lawhead está limpio, pero se me ocurrió la siguiente solución: usar matrices enmascaradas Numpy.
No sé si es mejor, pero eso era lo que sabía antes (hace 3 años ...).
Originalmente, creé un área de datos válida dentro del ráster original (p. Ej., La extensión del ráster de salida donde era igual), pero me gustó la idea de hacer que el ráster también sea más pequeño (p. Ej., Crop_to_cutline), así que adopté world2Pixel de Joel Lawhead . Aquí está mi propia solución:

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)

para obtener una descripción completa de class MaskRaster y sus métodos, consulte github de mi proyecto .

Usando este código, aún necesitarás usar GDAL. Sin embargo, el plan es usar en el futuro Python puro donde pueda, porque la audiencia prevista de mi software tiene dificultades con demasiadas dependencias (uso Debian para desarrollar el software y los clientes usan Windows 7 ...).

    
respondido por el Oz123 01.07.2012 - 12:29

Lea otras preguntas en las etiquetas