jueves, 5 de septiembre de 2019

Cómo hacer clipping de un ráster por rasgos individuales de un shapefile tipo polígono con GDAL/OGR en Python

En este post se incluye un procedimiento para hacer clipping de un ráster por rasgos individuales de un shapefile tipo polígono mediante los módulos Python de GDAL/OGR. La rutina introduce una proyección en las capas vectoriales de corte (que se eliminan al final de la rutina) e incluye la asignación de los nodata a -999 en las capas ráster (cuya proyección es tomada automáticamente del ráster base) dentro de la función gdal.Warp. El código completo es el siguiente:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
from osgeo import gdal, ogr, osr
import os

shape_file = "/home/zeito/pyqgis_data/boxes.shp"

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()
nr_features = layer.GetFeatureCount()

spatialRef = layer.GetSpatialRef()
wkt = spatialRef.ExportToWkt()
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)

spatialRef = osr.SpatialReference()
spatialRef.ImportFromWkt(wkt)

spatialRef.MorphToESRI()

shp_paths = []

for feat_id in range(nr_features):
    feature = layer.GetFeature(feat_id)
    id_no = feature.GetField("FID")
    outShapefile = os.path.join("/home/zeito/pyqgis_data/", "feature" + str(feat_id)+".shp")
    shp_paths.append(outShapefile)
    outPrj = os.path.join("/home/zeito/pyqgis_data/", "feature" + str(feat_id)+".prj")

    outDriver = ogr.GetDriverByName("ESRI Shapefile")

    # Remove output shapefile if it already exists
    if os.path.exists(outShapefile):
        outDriver.DeleteDataSource(outShapefile)
        
    # Create the output shapefile
    outDataSource = outDriver.CreateDataSource(outShapefile)
    outLayer = outDataSource.CreateLayer("feature" + str(feat_id), geom_type=ogr.wkbPolygon)
    
    file = open(outPrj, 'w')
    file.write(spatialRef.ExportToWkt())
    file.close()

    # Add an ID field
    idField = ogr.FieldDefn("id", ogr.OFTInteger)
    outLayer.CreateField(idField)
    
    # Create the feature and set values
    featureDefn = outLayer.GetLayerDefn()
    feat = ogr.Feature(featureDefn)
    geom = feature.GetGeometryRef()

    feat.SetGeometry(geom)
    feat.SetField("id", id_no)
    outLayer.CreateFeature(feat)
    feat = None
    
    # Save and close DataSource
    outDataSource = None

    OutTile = gdal.Warp("/home/zeito/pyqgis_data/cut" + str(feat_id) + ".tif", 
                        "/home/zeito/pyqgis_data/utah_demUTM12.tif", 
                        cutlineDSName="/home/zeito/pyqgis_data/cut" + str(feat_id) + ".shp",
                        cropToCutline=True,
                        dstNodata = -999)

dataSource = None

for path in shp_paths:
    outDriver.DeleteDataSource(path)

La situación inicial luce como sigue:


Después de ejecutado el script en la Python Console de QGIS 3.8, así se observan 2 de los ráster individuales al ser cargados en la Map View. El tercero se deja sin visualizar para corroborar que la individualización de rasgos fue ejecutada correctamente.



Publicado por José Guerrero en 23:37

No hay comentarios:

Publicar un comentario

Suscribirse a: Enviar comentarios (Atom)