Introduzione

Mi è capitato di dover riproiettare dei layer in QGIS, ho spesso usufruito della comoda riproiezione al volo oppure ho salvato il dataset in un diverso CRS. La situazione che mi ha portato alla scrittura dello script che andrò a dettagliare molto velocemente nel seguito è la seguente: un cliente aveva un progetto con centinaia di layer, sia vettoriali che raster, e aveva la necessità di riportare i soli layer vettoriali in un unico sistema di riferimento (alcuni lo erano già, altri no). Non si sapeva però quali fossero i layer che avevano un sistema di riferimento diverso da quello ricercato. I layer da controllare erano parecchi, motivo per cui ho prefereito usare alcuni semplici comandi di python da lanciare nella console che QGIS mette a disposizione, per automatizzare il processo di ricerca dei layer vettoriali con sistema di riferimento diverso da quello ricercato e riproiezione in una cartella da me specificata. Lo script finale che ho fatto girare è molto simile a quello riportato di seguito:

canvas = qgis.utils.iface.mapCanvas()
allLayers = canvas.layers()

pathToProject=os.getcwd()

if not os.path.exists('Reprojected'):
    os.makedirs('Reprojected')

outputCRSSstr = 'EPSG:3003'
outputCRS = QgsCoordinateReferenceSystem(outputCRSStr)

for i in allLayers: 
    if i.type() != QgsMapLayer.VectorLayer :

        print(i.name() + " skipped as it is not a vector layer")

    if ( (i.type() == QgsMapLayer.VectorLayer) and i.crs() != outputCRS ):

        print (i.name())

        outputLayerPathReproj = os.path.join(pathToProject, "Reprojected", i.name()+".shp")

        parameter = {'INPUT': i, 'TARGET_CRS': outputCRSStr, 'OUTPUT': '{}'.format(outputLayerPathReproj)}
        result = processing.run('native:reprojectlayer', parameter)

        vlayer = QgsVectorLayer(outputLayerPathReproj, "reproj"+ i.name(), "ogr")
        QgsProject.instance().addMapLayer(vlayer)

Non dettaglierò tutti i comandi usati comunemente in PyQGIS in quanto la documentazione ufficiale è completa e ben gestita, cercherò di spiegare come ho unito tutti i pezzi e come fare a ottenere determinati risultati; in sintesi: si può usare la documentazione ufficiale come un bignami per comporre lo script adatto alle proprie esigenze. In questo breve articolo mostrerò un'applicazione particolare che unisce diversi mattoncini presi dalla documentazione e che mi è servita per risolvere alcuni problemi durante una consulenza.

Spiegazione dello script

Le prime due righe ci permettono di accedere al map canvas e in particolare di accedere ai layer attivi nel progetto.

canvas = qgis.utils.iface.mapCanvas()
allLayers = canvas.layers()

Ho inizializzato due variabili per chiarezza del lettore, avrei anche potuto dichiare solo:

allLayers = qgis.utils.iface.mapCanvas().layers()

Il comando os.getcwd() restituisce la working directory di un processo come si può vedere dalla documentazione del modulo OS.

pathToProject=os.getcwd()

Mi servo del modulo OS per capire in che cartella sto lavorando per poter creare una nuova cartella, che chiamo Reprojected; prima di creare questa cartella nella stessa posizione da cui si fa partire lo script (la posizione del progetto .qgs) controllo che non esista già con un if statement, se non esiste la creo con os.makedirs:

if not os.path.exists('Reprojected'):
    os.makedirs('Reprojected')

A questo punto dichiaro il CRS in cui si vogliono ottenere i dati e lo definisco attraverso la classe QgsCoordinateReferenceSystem:

outputCRSSstr = 'EPSG:3003'
outputCRS = QgsCoordinateReferenceSystem(outputCRSStr)

Inizio ad analizzare ogni layer contenuto nella variabile allLayers attraverso un ciclo for che itera su tutti i suoi elementi. Dal momento che sono interessato ai soli layer vettoriali se il type() di uno degli elementi analizzati non è vettoriale lo salto nella riproiezione:

 for i in allLayers: 
   if i.type() != QgsMapLayer.VectorLayer :
       print(i.name() + " skipped as it is not a vector layer")

Se invece il layer in considerezione i è di tipo vettoriale e se il suo CRS è diverso dal CRS target proseguo con l'analisi

   if ( (i.type() == QgsMapLayer.VectorLayer) and i.crs() != outputCRS ):
       print (i.name())

Quindi:

• Definisco il percorso e il nome del layer riproiettato che si chiamerà con il nome del layer originale, e sarà contenuto nella cartella Reprojected che abbiamo creato in precedenza;
• Indico i parametri da dare in pasto all'algoritmo di processing che utilizzo per la riproiezione e lancio l'algoritmo stesso richiamando il suo identificativo (reprojectLayer) e indicandone i parametri;
• Infine carico i vettori riproiettati nel progett QGIS aperto con il nome uguale a quello dei layer originali ma con anteposta la stringa reproj;

        outputLayerPathReproj = os.path.join(pathToProject, "Reprojected", i.name()+".shp")

        parameter = {'INPUT': i, 'TARGET_CRS': outputCRSStr, 'OUTPUT': '{}'.format(outputLayerPathReproj)}
        result = processing.run('native:reprojectlayer', parameter)

        vlayer = QgsVectorLayer(outputLayerPathReproj, "reproj"+ i.name(), "ogr")
        QgsProject.instance().addMapLayer(vlayer)

Conclusione

A questo punto se il progetto QGIS da cui ho fatto partire lo script conteneva layer vettoriali con CRS diverso da quello target, dovrei ritrovare nello stesso progetto i layer riproiettati. Ovviamente lo script è costumizzabile e migliorabile, deve essere cambiato secondo le esigenze (ad esempio a seconda del CRS che si vuole ottenere) ma può essere uno spunto per gestire i propri dati in maniera ordinata e riproducibile da tutti.