Ortofoton som Cloud Optimized GeoTIFF (COG) från Lantmäteriet

[lmdaniel]

Numera finns ortofoton som avgiftsfria data från Lantmäteriet. Dessa är i sin tur tillgängliga som Cloud Optimized GeoTIFF (COG). Det ger möjligheten att strömma dessa i exempelvis QGIS.

Nedan finns ett pythonskript (QGIS) som inom kartfönstrets utbredningsområde (bbox) gör en sökning mot Lantmäteriets api och lägger till de senaste ortofotona inom området som strömmade bilder (vsicurl). Skriptet kan ju säkert optimeras mer.

import requests
import json
from qgis.core import QgsVectorLayer, QgsProject, QgsFeature, QgsGeometry, QgsPointXY, QgsField, QgsCoordinateTransform, QgsCoordinateReferenceSystem, QgsRasterLayer
from PyQt5.QtCore import QVariant
from qgis.utils import iface
import os
from requests.auth import HTTPBasicAuth
from pyproj import Transformer
from datetime import datetime

# Autentisering
session = requests.Session()
session.auth = HTTPBasicAuth("användarnamn", "lösenord")

# Tidsstämpel
timestamp = datetime.now().strftime("%Y%m%d%H%M%S")

# Hämta bbox i SWEREF 99 TM och konvertera till WGS 84
extent = iface.mapCanvas().extent()
bbox_sweref99 = [extent.xMinimum(), extent.yMinimum(), extent.xMaximum(), extent.yMaximum()]
transformer_to_wgs84 = Transformer.from_crs("EPSG:3006", "EPSG:4326")
bbox_wgs84 = [transformer_to_wgs84.transform(bbox_sweref99[1], bbox_sweref99[0]), transformer_to_wgs84.transform(bbox_sweref99[3], bbox_sweref99[2])]
bbox_wgs84 = [bbox_wgs84[0][1], bbox_wgs84[0][0], bbox_wgs84[1][1], bbox_wgs84[1][0]]

# Hämta data från API
api_url = "https://api.lantmateriet.se/stac-bild/v1/search"
response = session.post(api_url, headers={"Content-Type": "application/json"}, data=json.dumps({"bbox": bbox_wgs84, "limit": 100}))
features = response.json().get("features", [])

# Skapa lager i QGIS
def create_layer(name, features, crs="EPSG:3006"):
    layer = QgsVectorLayer(f"Polygon?crs={crs}", name, "memory")
    provider = layer.dataProvider()
    if features:
        first_feature = features[0]
        for key in first_feature['properties'].keys():
            provider.addAttributes([QgsField(key, QVariant.String)])
        for asset_key in first_feature['assets'].keys():
            provider.addAttributes([QgsField(f"asset_{asset_key}", QVariant.String)])
        layer.updateFields()
        transformer = QgsCoordinateTransform(QgsCoordinateReferenceSystem("EPSG:4326"), QgsCoordinateReferenceSystem(crs), QgsProject.instance())
        for feature in features:
            qgs_feature = QgsFeature()
            attributes = [str(feature['properties'].get(key, "")) for key in first_feature['properties'].keys()]
            for asset_key in first_feature['assets'].keys():
                attributes.append(feature['assets'][asset_key].get('href', ''))
            qgs_feature.setAttributes(attributes)
            geom = feature.get("geometry")
            if geom and geom.get("type") == "Polygon":
                points = [transformer.transform(QgsPointXY(coord[0], coord[1])) for coord in geom.get("coordinates")[0]]
                qgs_feature.setGeometry(QgsGeometry.fromPolygonXY([points]))
                provider.addFeatures([qgs_feature])
    QgsProject.instance().addMapLayer(layer)

create_layer("Utbredning ortofoton, alla år", features)

# Skapa lager med senaste året
latest_year = max(int(f['properties'].get('flygar', 0)) for f in features)
latest_features = [f for f in features if int(f['properties'].get('flygar', 0)) == latest_year]
create_layer(f"Utbredning ortofoton, {latest_year}", latest_features)

# Funktion för att skapa ett rasterlager från en URL
def create_raster_layer(url):
    vsicurl_url = f"/vsicurl/{url}"
    layer_name = os.path.basename(url).replace('.tif', '')
    raster_layer = QgsRasterLayer(vsicurl_url, layer_name, "gdal")
    if raster_layer.isValid():
        QgsProject.instance().addMapLayer(raster_layer)
        print(f"Lagret {layer_name} har lagts till i projektet.")
    else:
        print(f"Lagret {layer_name} kunde inte läsas in!")

# Skapa rasterlager för varje asset-länk (exkludera thumbnails)
for feature in latest_features:
    for asset_key, asset_info in feature['assets'].items():
        href = asset_info.get('href', '')
        if href and not href.endswith('thumbnail.jpg'):
            create_raster_layer(href)

Har du några egna erfarenheter att jobba med Cloud Optimized GeoTIFF?