QGis e PythonPath

Ultimamente ho incasinato la Debian box e mi compare questo errore aprendo QGis...dopo funziona ma non disponibili i plugin 

la soluzione e' settare in modo esplicito la PythoPath

PYTHONPATH=/usr/lib/python3/dist-packages:$PYTHONPATH 

qgis

Enmap Toolbox e Debian 13

Niente da fare ...Enmap Toolbox plugin di QGis non funziona con Debian 13...ci ho provato in tutti i modi (il problema principale ad una libreria di Qt5) ed alla fine sono tornato a montare una macchina virtuale di Debian 12 che installa QGis 3.22 LTS  

 

  

python3 -m venv virtualenv --system-site-packages

source virtualenv/bin/activate 

pip install enpt_enmapboxapp

qgis 

Anomalie in NDVI con algoritmo Continuous Change Detection (CCDC)

Attenzione: questo plugin non funziona su Debian 13 per delle dipendenze non risolte. Funziona in modo corretto su  Debian 12 con QGis 3.22 LTS

Attenzione 2 : la funzione CCDC e' presente su Google Earth Engine solo sui profili a pagamento...per questo motivo ho usato la versione di QGis...Aggiornamento: a quanto sembra non e' necessario il profilo a pagamento ma che lo script sia collegato ad un progetto di Google Cloud ed autorizzato (un progetto autorizzato si riconosce dall'icona in alto a destra)

Progetto autorizzato

Progetto non autorizzato

 

 Questo plugin usa l'algorimo Continous Change Detection che si basa su una funzione armonica per evidenziare brusche variazioni del parametro analizzato

Questo e' lo script per Google Earth Engine

// 1. DEFINIZIONE ROI

var areaOfInterest = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[10.896754172154326, 43.02122523564496],
          [10.896754172154326, 43.002271538414895],
          [10.935463812657256, 43.002271538414895],
          [10.935463812657256, 43.02122523564496]]], null, false);

Map.centerObject(areaOfInterest, 15);

// 2. PREPARAZIONE DATI
var prepareData = function(image) {
  var optical = image.multiply(0.0000275).add(-0.2);
  var ndvi = optical.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI');
  return image.addBands(ndvi).select('NDVI');
};

var collection = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2")
  .merge(ee.ImageCollection("LANDSAT/LT05/C02/T1_L2"))
  .merge(ee.ImageCollection("LANDSAT/LE07/C02/T1_L2"))
  .filterBounds(areaOfInterest)
  .map(prepareData)
  .sort('system:time_start');

// 3. ESECUZIONE CCDC CON GESTIONE ERRORI
// Proviamo a chiamare l'algoritmo nativo in modo più "protetto"
var ccdc;
try {
  ccdc = ee.Algorithms.TemporalSegmentation.ccdc({
    collection: collection,
    breakThreshold: 1.0,
    minObservationsForSegment: 6
  });
  print('CCDC eseguito con successo:', ccdc);
} catch (e) {
  print('Errore critico: Il tuo account non ha accesso a CCDC tramite API diretta.', e);
  print('Soluzione: Assicurati che lo script sia collegato a un Cloud Project (in alto a destra).');
}

// 4. SE CCDC HA FUNZIONATO, VISUALIZZA
if (ccdc) {
  var magnitude = ccdc.select('NDVI_magnitude').arrayGet([0]);
  var tBreak = ccdc.select('tBreak').arrayGet([0]);

  var magViz = {min: -0.5, max: 0.5, palette: ['#d73027', '#f4a582', '#f7f7f7', '#d9f0d3', '#1a9850']};
  var dateViz = {min: 1984, max: 2024, palette: ['#0000ff', '#00ffff', '#ffff00', '#ff0000']};

  Map.addLayer(magnitude.clip(areaOfInterest), magViz, 'Magnitudo Cambiamento');
  Map.addLayer(tBreak.clip(areaOfInterest), dateViz, 'Anno del Cambiamento');
}

Preparazione dati per segmentazione rete neurale

 Ho trovato in rete un dataset di foto da drone con gia' classificate le fratture. Si tratta di dati con licenza permissiva 

 

 

I dati derivano da questo articolo 

https://www.researchgate.net/publication/357994422_A_new_subsampling_methodology_to_optimize_the_characterization_of_two-dimensional_bedrock_fracture_networks

@article{article,
author = {Ovaskainen, Nikolas and Nordbäck, Nicklas and Skyttä, Pietari and Engström, Jon},
year = {2022},
month = {01},
pages = {104528},
title = {A new subsampling methodology to optimize the characterization of two-dimensional bedrock fracture networks},
volume = {155},
journal = {Journal of Structural Geology},
doi = {10.1016/j.jsg.2022.104528}
}

e sono disponibili come ortofoto RGB da drone a questo link https://zenodo.org/records/4719627  (16 Gb) mentre le fratture sono tematizzate in un diversi file gpkg a questo link https://www.kaggle.com/datasets/nasurz/getaberget-fracture-trace-dataset/data (4Mb)

il problema e' che in questo formato le immagini non sono adatte al processamento per la rete neurale. L'idea e' quindi di dividere le ortofoto in tasselli da 15 metri e creare una immagine maschera con le sole fratturazioni. Tale compito e' stato eseguito con uno script Python in QGis

Nello script si devono impostare le coordinate geografiche della finestra  ed il folder dove saranno salvate le varie tiles. In un primo passaggio si disattiva dal progetto QGis il tema lineare delle fratture per avere solo i tasselli delle ortofoto e dopo si disattiva il tematismo raster per creare le maschere

 ------------------------------------------------------------------------------------------

from qgis.core import (
    QgsProject, QgsApplication, QgsLayout, QgsPrintLayout, QgsLayoutItemMap,
    QgsUnitTypes, QgsLayoutSize, QgsLayoutPoint, QgsLayoutExporter, QgsRectangle
)
import os



# Path to your QGIS project
project_path = "C:\\Users\\l.innocenti\\Desktop\\geta.qgz"

# Output folder
output_dir = "C:\\Users\\l.innocenti\\Desktop\\geta5\\mask\\"
os.makedirs(output_dir, exist_ok=True)

# Load project
project = QgsProject.instance()
project.read(project_path)


from qgis.core import QgsRectangle

# Compute project extent from all layers
extent = None
for layer in project.mapLayers().values():
    if not layer.isValid():
        continue
    layer_extent = layer.extent()
    if extent is None:
        extent = QgsRectangle(layer_extent)
    else:
        extent.combineExtentWith(layer_extent)

if extent is None:
    raise ValueError("No valid layers found in the project.")

xmin, ymin = 108914, 6720011
xmax, ymax = 108969, 6720066

# Get the full extent in project CRS
#extent = project.boundingBox()
#xmin, ymin, xmax, ymax = extent.xMinimum(), extent.yMinimum(), extent.xMaximum(), extent.yMaximum()

# Tile size in meters
tile_size = 15

# Iterate over the grid
x = xmin
tile_id = 0
while x < xmax:
    y = ymin
    while y < ymax:
        # Tile extent
        tile_extent = QgsRectangle(x, y, x + tile_size, y + tile_size)

        # Create a print layout
        layout = QgsPrintLayout(project)
        layout.initializeDefaults()
        layout.setName(f"tile_{tile_id}")
        pc = layout.pageCollection().pages()[0]

        # Map item
        map_item = QgsLayoutItemMap(layout)
        map_item.setRect(20, 20, 200, 200)
        map_item.setExtent(tile_extent)
        map_item.setBackgroundColor(QColor(255, 255, 255, 0))
        map_item.attemptResize(QgsLayoutSize(200, 200, QgsUnitTypes.LayoutMillimeters))
        visible_layers = [lyr for lyr in iface.mapCanvas().layers()]
        map_item.setLayers(visible_layers)

        layout.addLayoutItem(map_item)
        layout.addLayoutItem(map_item)
        map_item.attemptMove(QgsLayoutPoint(0, 0, QgsUnitTypes.LayoutMillimeters))

        # Export PNG
        exporter = QgsLayoutExporter(layout)
        out_path = os.path.join(output_dir, f"tile_{tile_id}.png")
        exporter.exportToImage(out_path, QgsLayoutExporter.ImageExportSettings())

        print(f"Saved {out_path}")
        tile_id += 1

        y += tile_size
    x += tile_size
  ------------------------------------------------------------------------------------------

 alla fine si hanno coppie di questo tipo (le immagini sono state rifilate con imagemagick per togliere lo spazio bianco inserito da QGis)

 

 

 

 per "aggiustare" il background con colore nero e il target come grigio come 127 

import cv2
import numpy as np
import sys

# --- Check for input filename ---
if len(sys.argv) < 2:
    print("Usage: python script.py <image_filename>")
    sys.exit(1)

filename = sys.argv[1]

# Read image
img = cv2.imread(filename)

if img is None:
    raise FileNotFoundError(f"Cannot read file: {filename}")

# Define white color in BGR
white = np.array([255, 255, 255], dtype=np.uint8)

# Create mask for white pixels
mask_white = np.all(img == white, axis=-1)

# Create output image
output = np.full_like(img, (128, 128, 128))  # Set all to gray
output[mask_white] = (0, 0, 0)  # Set white pixels to black

# Save with same filename (overwrite original)
cv2.imwrite(filename, output)

print(f"Image processed and saved: {filename}")

 

 

Enmap ToolBox su Linux

Per installare Enmap ToolBox su Debian Linux si devono aggiungere alcune dipendenze. Alcune sono presenti gia' in apt ma enpt_enmapboxapp deve essere scaricata con pip. A questo punto e' necessario creare un venv e si deve lanciare qgis dall'interno del venv

sudo apt install python3-h5py python3-pyqt5.qtopengl python3-netcdf4

python3 -m venv .venv --system-site-packages

pip3 install enpt_enmapboxapp

Plugin OrfeoToolbox in QGis Windows

 Per utilizzare Orfeo Toolbox direttamente dentro QGis si spacchetta lo zip di Orfeo Toolbox in una directory (con il nome senza caratteri speciali). Poi si seleziona strumenti di processing/Options /Programmi e si popolano i campi 

Cartella OTB : C:/OTB-7.4.1-Win64/OTB-7.4.1-Win64

Cartella delle applicazioni OTB : C:/OTB-7.4.1-Win64/OTB-7.4.1-Win64/bin;C:/OTB-7.4.1-Win64/OTB-7.4.1-Win64/lib/otb/applications

Pansharpening Landsat con Qgis

 

Si caricano i raster delle bande 2,3,4 (BGR) e la banda 8 pancromatica

Si apre il menu Raster/Miscellanea/Crea Raster Virtuale selezionando come input le bande ottiche e spuntando il flag Place each input file in a separate band

A questo punto dal toolbox si seleziona pansharpening indicando come Dataset spettrale il raster virtuale creato in precedenza e come dataset pancromatico la banda 8

Intersezione superficie topografica con piani con qgSurf

Un plugin di QGis che non conoscevo ma che puo' tornare utile per verificare la correttezza delle faglie disegnate su una carta geologica.

Questo plugin, presi un modello digitale del terreno e un piano, calcola la linea di intersezione tra le due superfici...una regola della V tecnologica

Per tentativi puo' essere anche stimata la pendenza del piano di faglia ove sia riesca a mappare il contatto in campagna

Stereoscopio interattivo con QGis2ThreeJS

Circa 3 anni fa avevo provato il plugin di QGis Qgis2ThreeJS per fare DEM stereoscopici da usare Google CardBoard ma poi avevo lasciato li' la cosa perche' non riuscivo ad avere l'interazione per girare il modello. Adesso ho trovato la soluzione ed e' molto piu' banale di quanto potessi immaginare...un mouse bluetooth. Basta accoppiare il mouse al telefono, inserire il telefono in un CardBoard o DayDream e si puo' visualizzare il modello da piu' punti di vista

Per fare il progetto QGis si deve utilizzare il plugin nel ramo sperimentale

https://github.com/minorua/Qgis2threejs/tree/exp_stereo

e si deve selezione il template 3DViewer.html

Si seleziona poi il modo di controllare il modello

e questo il risultato finale (ovviamente da inserire nel visore)

QGIS2ThreeJs

QGis2ThreeJS e' un comodo plugin di QGis per creare dei modelli digitali dal terreno che possono essere visualizzati in modo interattivo su Web mediante la libreria ThreeJS

Per la prova ho preso i dati DEM ripresi da Aster (http://earthexplorer.usgs.gov/)

Non avendo delle immagini georiferite da spalmare sul DEM ho ripreso una Geotiff del rilievo a falsi colori da questo indirizzo http://geodati.fmach.it/gfoss_geodata/SRTM-Italy/

Una volta caricati i dati come raster su QGis si apre il plugin QGis2ThreeJS  scegliendo il tema del DEM ed il tema dell'immagine

deve essere selezionata anche l'esagerazione verticale (a tentativi ho trovato il valore di 0.02)

Una volta scelto il nome del file del progetto html si clicca Run e si apre in automatico il browser e vengono caricati i dati generati dal plugin

Problemi con QGis in Debian

Al momento di usare QGis su Debian mi e' comparso il seguente messaggio

La cosa curiosa e' che avevo di recente usato QGis senza errori di sorta. Ho verificato che il malfunzionamento deriva dall'aggiornamento di alcune librerie di appoggio di QGis e del mancato aggiornamento di QGis stesso (il repository ufficiale non e' sincronizzato)
La soluzione e' quella di impiegare come repository direttamente qgis.org

per prima cosa si rimuove la versione di default di QGis

apt-get remove qgis python-qgis qgis-plugin-grass

si aggiunge al file /etc/apt/sources.list il repository di QGis

deb http://qgis.org/debian jessie main

si aggiungono le chiavi
gpg --recv-key 47765B75 gpg --export --armor 47765B75 | sudo apt-key add -
e si reinstalla QGis dalla nuova sorgente

sudo apt-get update sudo apt-get install qgis python-qgis qgis-plugin-grass

come risultato sono passato da una versione 1.7 non funzionante (o meglio con i plugin Python non funzionanti) ad una fiammante versione 2.0

Cambio sistema di coordinate geografiche

Ogni tanto mi ricordo di essere un geologo con un po' di esperienza in telerilevamento.
Per un studio che sto facendo nei ritagli di tempo sto usando delle immagini Landsat 7 e non avendo una licenza Envi e ArcGis mi sto adattando ad usare QGis

Caricando i dati Landsat Toscana-Emilia (TIF georiferite) direttamente in QGis queste vengono inserite nel sistema di riferimento EPSG:32632 (WGS84 UTM Zona 32 N)

LE71920292001046EDC01

Gli shapefile che invece mi servono sono tutti nel sistema di riferimento EPSG 3003 (Monte Mario Italy Zona 1)

Invece di convertire le immagini e' possibile traslare il sistema di riferimento degli shapefile usando i tools GDAL in particolare con il comando

dove
wgs.shp e' il nome dello shape risultato dell'elaborazione
boaga.shp e' il nome del file shape da convertire
-s_srs e' il sistema di riferimento di partenza
-t_srs e' il sistema di riferimento finale

ogr2ogr wgs.shp boaga.shp -s_srs EPSG:3003 -t_srs EPSG:32632