Embedded Files
Per scaricare il file a pezzi seguire la procedura
Per scaricare il file a pezzi seguire la procedura
apri Qgis
carica una mappa di base
apri in Plugins consol Python - mostra editor
copia nell'editor il seguente script (file txt allegato di seguito)
premi Run
fa una finestra sulla mappa caricata tenedo premito il tasto del mouse e rilasciandolo nel vertice opposto.
# -*- coding: utf-8 -*-
"""
draw_bbox_clip_dtm.py
=====================
Da incollare nella Console Python di QGIS.
Attiva uno strumento sul map canvas: si traccia un rettangolo trascinando il
mouse, e al rilascio viene scaricato SOLO il DTM contenuto nella bbox disegnata,
via GDAL /vsicurl/. Il ritaglio viene poi caricato come layer nel progetto.
Funziona perche' /vsicurl/ legge via HTTP Range solo i tile che coprono l'area:
da un DTM nazionale di decine di GB si estraggono pochi MB.
Configura le due variabili qui sotto e lancia. Poi disegna il rettangolo.
"""
from qgis.core import (
QgsRectangle, QgsPointXY, QgsWkbTypes, QgsProject, QgsRasterLayer,
QgsCoordinateReferenceSystem, QgsCoordinateTransform,
)
from qgis.gui import QgsMapTool, QgsRubberBand
from qgis.PyQt.QtGui import QColor
from qgis.utils import iface
from osgeo import gdal
import os
# ============================ CONFIGURAZIONE ============================
# URL del DTM remoto (GeoTIFF/COG raggiungibile via HTTP/HTTPS).
# DTM 5 m d'Italia (HR-DTM-5m, Zenodo).
DTM_URL = "https://zenodo.org/records/18872933/files/HRDTM5m"
# Cartella di output per il ritaglio.
OUT_DIR = os.path.expanduser("~/dtm_clip")
# =======================================================================
# Suggerimenti di rete per /vsicurl/ (meno richieste HTTP inutili)
gdal.SetConfigOption("GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR")
gdal.SetConfigOption("GDAL_HTTP_MULTIPLEX", "YES")
gdal.SetConfigOption("VSI_CACHE", "TRUE")
def _raster_crs_and_type(vsi_path):
"""Apre il raster remoto (solo header) e ne restituisce CRS e dtype."""
ds = gdal.Open(vsi_path)
if ds is None:
raise RuntimeError("Impossibile aprire il raster remoto: " + vsi_path)
wkt = ds.GetProjection()
crs = QgsCoordinateReferenceSystem.fromWkt(wkt)
band = ds.GetRasterBand(1)
is_float = gdal.GetDataTypeName(band.DataType).lower().startswith("float")
ds = None
return crs, is_float
def clip_dtm_to_bbox(rect_canvas, canvas_crs):
"""Ritaglia il DTM remoto sulla bbox (nel CRS del canvas) e carica il layer."""
vsi_path = "/vsicurl/" + DTM_URL
raster_crs, is_float = _raster_crs_and_type(vsi_path)
# Riproietta la bbox dal CRS del canvas a quello del raster, se diversi
if canvas_crs != raster_crs:
xform = QgsCoordinateTransform(canvas_crs, raster_crs, QgsProject.instance())
rect = xform.transformBoundingBox(rect_canvas)
else:
rect = rect_canvas
os.makedirs(OUT_DIR, exist_ok=True)
out_path = os.path.join(OUT_DIR, "dtm_aoi.tif")
# PREDICTOR=3 per i dati float (DTM), 2 per gli interi
predictor = "3" if is_float else "2"
# projWin vuole [ulx, uly, lrx, lry] = [xmin, ymax, xmax, ymin]
print("==> Ritaglio in corso (scarico solo i tile dell'area)...")
gdal.Translate(
out_path,
vsi_path,
projWin=[rect.xMinimum(), rect.yMaximum(),
rect.xMaximum(), rect.yMinimum()],
creationOptions=["COMPRESS=DEFLATE", "PREDICTOR=" + predictor, "TILED=YES"],
)
layer = QgsRasterLayer(out_path, "DTM AOI")
if not layer.isValid():
raise RuntimeError("Il ritaglio non e' un layer valido: " + out_path)
QgsProject.instance().addMapLayer(layer)
print("==> Fatto: " + out_path)
size_mb = os.path.getsize(out_path) / 1e6
print(" Dimensione scaricata: {:.1f} MB".format(size_mb))
class RectangleClipTool(QgsMapTool):
"""Strumento di disegno: trascina per definire un rettangolo, rilascia per ritagliare."""
def __init__(self, canvas):
super().__init__(canvas)
self.canvas = canvas
self.rubber = QgsRubberBand(canvas, QgsWkbTypes.PolygonGeometry)
self.rubber.setColor(QColor(193, 36, 24, 60)) # rosso Planetek, riempimento
self.rubber.setStrokeColor(QColor(193, 36, 24)) # bordo pieno
self.rubber.setWidth(2)
self.start_point = None
self.dragging = False
def _draw_rect(self, p1, p2):
rect = QgsRectangle(p1, p2)
self.rubber.reset(QgsWkbTypes.PolygonGeometry)
corners = [
QgsPointXY(rect.xMinimum(), rect.yMinimum()),
QgsPointXY(rect.xMaximum(), rect.yMinimum()),
QgsPointXY(rect.xMaximum(), rect.yMaximum()),
QgsPointXY(rect.xMinimum(), rect.yMaximum()),
]
for i, pt in enumerate(corners):
self.rubber.addPoint(pt, i == len(corners) - 1)
self.rubber.show()
def canvasPressEvent(self, event):
self.start_point = self.toMapCoordinates(event.pos())
self.dragging = True
def canvasMoveEvent(self, event):
if self.dragging and self.start_point is not None:
self._draw_rect(self.start_point, self.toMapCoordinates(event.pos()))
def canvasReleaseEvent(self, event):
if not self.dragging:
return
self.dragging = False
end_point = self.toMapCoordinates(event.pos())
rect = QgsRectangle(self.start_point, end_point)
rect.normalize()
self.rubber.reset(QgsWkbTypes.PolygonGeometry)
canvas_crs = self.canvas.mapSettings().destinationCrs()
self.canvas.unsetMapTool(self)
if rect.width() == 0 or rect.height() == 0:
print("Rettangolo nullo: riprova trascinando il mouse.")
return
try:
clip_dtm_to_bbox(rect, canvas_crs)
except Exception as e:
print("Errore durante il ritaglio: " + str(e))
# Attiva lo strumento. NB: va tenuto un riferimento globale (_bbox_tool),
# altrimenti Python lo distrugge e lo strumento smette di funzionare.
_bbox_tool = RectangleClipTool(iface.mapCanvas())
iface.mapCanvas().setMapTool(_bbox_tool)
print("Strumento attivo: traccia un rettangolo sul canvas per ritagliare il DTM.")
Page updated
Google Sites
Report abuse