Buen día, Necesito ayuda para el código. Quiero realizar un script para un recorte masivo de polígonos de una capa vectorial sobre un raster. Dispongo de un shapefile con, por ejemplo 10 polígonos distintos, y quiero que este shapefile recorte mi capa raster y obtenga como resultado, 10 imágenes independientes que corresponda cada una de ellas a un polígono distinto. En la consola Python de Qgis consigo mediante el código:
from osgeo import gdal
import os
from pathlib import Path
basepath = Path().home().joinpath("OneDrive", "Escritorio", "Presentacion_Script")
rasterPath = basepath.joinpath("Ortos")
vectorPath = basepath.joinpath("Shapes", "Conjunto_1.shp")
#Guarda el recorte en una carpeta
dir_recorte = rasterPath.joinpath('Carpeta recortes')
if not os.path.exists(dir_recorte):
os.mkdir(dir_recorte)
#Recorte de la imagen
for lyr in rasterPath.glob("*.tif"):
output_path = dir_recorte.joinpath("ID_uno_{}".format(lyr.name))
parameters = {
'ALPHA_BAND': False,
'CROP_TO_CUTLINE': True,
'DATA_TYPE': 0,
'INPUT': lyr.as_posix(),
'KEEP_RESOLUTION': False,
'MASK': vectorPath.as_posix(),
'MULTITHREADING': False,
'OPTIONS': '',
'OUTPUT': output_path.as_posix(),
'SET_RESOLUTION': False,
'SOURCE_CRS': None,
'TARGET_CRS': None,
'X_RESOLUTION': None,
'Y_RESOLUTION': None
}
resultado = processing.run("gdal:cliprasterbymasklayer", parameters)
QgsProject.instance().addMapLayer(QgsRasterLayer(resultado['OUTPUT'], output_path.stem, 'gdal'))
Tras muchas pruebas intento realizar un bucle for para que recorte cada uno de los polígonos y obtener así las 10 imágenes que se necesitan.
Por ejemplo el shapefile orginal, con las 10 entidades, se ha dividido en un shapefile independiente cada una de las entidades para que el bucle, escoja dentro de la carpeta los archivos .shp y una iteración sobre estos archivos .shp, el script recorte los 10 shapefile, pero no consigo que el bucle funcione.
from osgeo import gdal
import os
from pathlib import Path
basepath = Path().home().joinpath("OneDrive", "Escritorio", "Doctorado", "Presentacion_Script")
rasterPath = basepath.joinpath("Ortos")
vectorPath = basepath.joinpath("Shapes", "Conjunto_1.shp")
ejemplo_dir = '/Users/Lenovo/OneDrive/Escritorio/Doctorado/Presentacion_script/Shapes/'
contenido = os.listdir(ejemplo_dir)
shape = []
for fichero in shape:
if os.path.isfile(os.path.join(ejemplo_dir, fichero)) and fichero.endswith('.shp'):
shape = []
shape.append(fichero)
with os.scandir(ejemplo_dir) as ficheros:
ficheros = [fichero.name for fichero in ficheros if fichero.is_file() and fichero.name.endswith('.shp')]
#Guarda el recorte en una carpeta
dir_recorte = rasterPath.joinpath('Carpeta recortes')
if not os.path.exists(dir_recorte):
os.mkdir(dir_recorte)
#Recorte de la imagen
for lyr in rasterPath.glob("*.tif"):
output_path = dir_recorte.joinpath("ID_uno_{}".format(lyr.name))
parameters = {
'ALPHA_BAND': False,
'CROP_TO_CUTLINE': True,
'DATA_TYPE': 0,
'INPUT': lyr.as_posix(),
'KEEP_RESOLUTION': False,
'MASK': vectorPath.as_posix(),
'MULTITHREADING': False,
'OPTIONS': '',
'OUTPUT': output_path.as_posix(),
'SET_RESOLUTION': False,
'SOURCE_CRS': None,
'TARGET_CRS': None,
'X_RESOLUTION': None,
'Y_RESOLUTION': None
}
resultado = processing.run("gdal:cliprasterbymasklayer", parameters)
QgsProject.instance().addMapLayer(QgsRasterLayer(resultado['OUTPUT'], output_path.stem, 'gdal'))
Tampoco he conseguido que el shapefile original sea recortado por cada uno de sus ID.
¿Me pueden ayudar para poder conseguir las 10 imágenes de una sola capa ráster?
Muchas gracias.
Juanam