Ritaglio delle mappe mediante shapefile

edigiacomo commented 3 years ago

Alcuni output richiedono il ritaglio della mappa mediante shapefile, ad esempio:

È una operazione che, a partire dal PNG generato da Magics, può essere effettuata con i seguenti comandi:

# Converto il PNG in GeoTiff perché gdalwarp non permette di indicare l'extent del file da riga di comando
# Ipotizzo che la mappa sia stata richiesta in EPSG:4326  e che `LONMIN LATMIN LATMIN LATMAX` 
# siano gli estremi del bbox.
# Se si usa il bbox espresso in WGS84 e la mappa è in una proiezione che non usa latlon (e.g. EPSG:3857)
# si devono riproiettare i punti.
gdal_translate -a_srs EPSG:4326 -a_ullr LONMIN LATMAX LONMAX LATMIN mappa.png mappa.tiff

# Genero un GeoTiff perché gdalwarp non può creare direttamente un PNG
gdalwarp -cutline emilia-romagna-4326.geojson mappa.tiff mappa-clip.tiff

# Creo il PNG
gdal_translate -of PNG mappa-clip.tiff mappa.png

Volendo, si può valutare di non implementare questa funzionalità in arkimaps ma di usare i suddetti comandi lato client.

Questa operazione ha senso (o quantomeno per come l'ho implementata io) solo con output georeferenziato (issue collegata #45).

spanezz commented 3 years ago

Se si può fare con chiamate a gdal_* si dovrebbe poter fare anche con python-gdal senza chiamare programmi esterni.

Abbiamo comunque come precondizione l'avere delle ricette che creano immagini georeferenziabili, senza bordini o legende all'interno. Immagino questa issue debba seguire #45

edigiacomo commented 2 years ago

Di seguito, un esempio di come è possibile fare il clip di un PNG usando uno shapefile senza usare file intermedi.

from osgeo import gdal, osr
# bounding box del PNG
bbox = [7, 43, 14, 46.5]
# dimensioni in pixel
size = [600, 422]
# Apro il file PNG e lo georeferenzio
ds_png = gdal.Open("input.png")
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds_png.SetProjection(srs.ExportToWkt())
ds_png.SetGeoTransform([
    bbox[0],
    (bbox[2]-bbox[0])/size[0],
    0,
    bbox[3],
    0,
    (bbox[1]-bbox[3])/size[1],
])
# Faccio il clip usando uno shapefile
gdal.Warp(
    "clip.png",
    ds_png,
    cutlineDSName="er.geojson",
    format="png"
)

Questa soluzione può essere applicata anche a #91 - che si era detto di considerare un caso particolare di questa issue.

edigiacomo commented 2 years ago

Dimenticavo i file con cui far girare l'esempio files.zip

spanezz commented 1 year ago

Ho iniziato il lavoro nel branch issue74. Apro alcuni elementi di discussione, intanto che continuo a sperimentare:

Proiezione dell'immagine PNG

Lo script di esempio ha hardcoded la georeferenziazione del PNG in input, perché i file PNG generati da Magics non sono georeferenziati.

Dai parametri di add_basemap possiamo tirar fuori subpage_map_projection e il bounding box dell'immagine, e sperare che siano sufficienti. Noto che lo script di esempio usa EPSG:4326 (WGS84), mentre emro_web e gli output tiled usano EPSG:3857 (web mercator), e intuisco che ci sarà da fare abbastanza attenzione a questo passaggio.

Caricare il geojson una volta sola

Caricare il geojson una volta per ogni immagine generata mi sembra inefficiente, e mi piacerebbe poterlo caricare, convertire nel sistema di coordinate dell'immagine da ritagliare, e riusare per tutte le immagini di un batch.

Non mi sono ancora chiari tutti i passaggi per farlo, ma ci si può arrivare.

Una volta ottenuta la geometria buona per il ritaglio, c'è da capire come passarla a gdal.Warp. In alternativa, cosa che forse può essere piú efficiente, si può sare gdal.Rasterize, una volta capito come.

Alcuni link:

Aggiungere "step" per il ritaglio

Gli step attuali sono tutti step di Magics, e l'immagine PNG vera e propria viene creata dopo che tutti gli step sono stati eseguiti da Magics. Come specifichiamo gli step di ritaglio?

Alcune opzioni:

Come step anche loro, che vengono riconosciuti come non-magics ed eseguiti alla fine. Bonus: il linter che segnala se gli step non-magics non vengono tutti dopo gli step magics
Come 'postprocess:' specificato come sequenza distinta di step
Solo nel flavour, e non nella ricetta

edigiacomo commented 1 year ago

Il seguente codice permette di fare il clipping del PNG "georeferenziabile" generato da arkimaps usando uno shapefile qualsiasi caricato in memoria e senza appoggiarsi a file temporanei salvati su disco.

Allego lo zip che contiene codice e file di input (PNG e GeoJSON) prototipo_issue74.zip, ma incollo anche qui sotto lo script nel caso in cui si volessero fare eventuali commenti a qualche riga.

Rispondo anche a qualche considerazione fatta:

Dai parametri di add_basemap possiamo tirar fuori subpage_map_projection e il bounding box dell'immagine, e sperare che siano sufficienti. Noto che lo script di esempio usa EPSG:4326 (WGS84), mentre emro_web e gli output tiled usano EPSG:3857 (web mercator), e intuisco che ci sarà da fare abbastanza attenzione a questo passaggio.

Nello script ho messo dei parametri di input che dovrebbero essere passati dalla ricetta. Nel caso della proiezione, si dovrà fare probabilmente un mapping tra alcune proiezioni "prosaiche" e l'EPSG corrispondente: sono comunque un numero limitato (26, si veda https://confluence.ecmwf.int/display/MAGP/Subpage+-+Projection) e non è detto che serva supportare tutte. Nel caso di 3857, abbiamo subpage_map_projection = EPSG:3857, quindi è facile :)

Caricare il geojson una volta per ogni immagine generata mi sembra inefficiente, e mi piacerebbe poterlo caricare, convertire nel sistema di coordinate dell'immagine da ritagliare, e riusare per tutte le immagini di un batch.

Nello script ho usato gdal.Rasterize (come suggerito da @spanezz) che prende un dataset vettoriale, che quindi non deve essere riletto tutte le volte da file. Non serve nemmeno convertirlo nel CRS dell'immagine da ritagliare - anche se forse potrebbe aver senso farlo per non delegarlo a ogni chiamata di gdal.Rasterize - ma lo terrei come ottimizzazione futura.

PNG input: input

PNG output: output

#!/usr/bin/python3
import osgeo
from osgeo import gdal, osr

# PARAMETRI FORNITI DALLA RICETTA

# bounding box del PNG in EPSG:4326 nel formato [LONMIN, LATMIN, LONMAX, LATMAX]
BBOX = [9.1951463, 43.7133281, 12.8283720, 45.1425750]
# dimensioni in pixel del PNF
SIZE = [1280, 705]
# epsg del PNG
EPSG = 3857

# FINE PARAMETRI FORNITI DALLA RICETTA

def convert_magics_bbox_to_epsg(bbox: tuple[float, float, float], epsg_out: int) -> tuple[float, float, float, float]:
    """Metodo per la conversione del bounding box di Magics (in EPSG:4326 e nella forma [LONMIN, LATMIN, LONMAX,
    LATMAX])."""
    srs_src = osr.SpatialReference()
    srs_src.ImportFromEPSG(4326)
    srs_dst = osr.SpatialReference()
    srs_dst.ImportFromEPSG(epsg_out)
    transform = osr.CoordinateTransformation(srs_src, srs_dst)
    # NOTA: in GDAL < 3 non viene rispettato l'ordine degli assi e viene sempre usato l'ordine lon,lat. Vedi
    # https://gdal.org/tutorials/osr_api_tut.html#crs-and-axis-order EPSG:4326 usa invece l'ordine lat,lon Se tutte le
    # operazioni successive con il bbox risultante sono delegate a GDAL (e.g. SetGeoTransform) non importa più l'ordine
    # in quanto è consistente con la versione di GDAL stessa.
    if osgeo.gdal_version[0] < 3:
        lr = transform.TransformPoint(bbox[0], bbox[1])
        ul = transform.TransformPoint(bbox[2], bbox[3])
    else:
        lr = transform.TransformPoint(bbox[1], bbox[0])
        ul = transform.TransformPoint(bbox[3], bbox[2])

    return [lr[0], lr[1], ul[0], ul[1]]

# Converto il bounding box della ricetta nella proiezione del raster
bbox = convert_magics_bbox_to_epsg(BBOX, EPSG)

# Apro il file PNG e lo georeferenzio
ds_png = gdal.Open("input.png")
srs = osr.SpatialReference()
srs.ImportFromEPSG(EPSG)
ds_png.SetProjection(srs.ExportToWkt())
ds_png.SetGeoTransform([
    bbox[0],
    (bbox[2]-bbox[0])/SIZE[0],
    0,
    bbox[3],
    0,
    (bbox[1]-bbox[3])/SIZE[1],
])
# Converto il PNG in un GeoTIFF in memoria usando il virtual FS di GDAL
# https://gdal.org/user/virtual_file_systems.html#vsimem-in-memory-files
# NOTA: si potrebbe anche usare il driver MEM e copiare a mano tutti i metadati e dati, ma non offre vantaggi se non
# quello di non dover fare l'unlink del file virtuale associato (vedi in fondo).
# NOTA: il file è visibile all'interno del processo, quindi in caso di multiprocesso non è un problema, ma nel caso di
# multithreading si deve generare un nome univoco (e.g. con threading.get_ident())
tif_path = "/vsimem/input.tif"
ds_tif = gdal.Translate(tif_path, ds_png, format="GTiff")
# Apro il file vettoriale
# NOTA: ogr.Open non va bene perché gdal.Rasterize vuole un tipo GDALDatasetShadow
# NOTA: Non serve riproiettare il vettoriale (che è in EPSG:32632, diversa dal raster)
shape = gdal.OpenEx("er.geojson", gdal.OF_VECTOR)
# Ritaglio il GeoTIFF, scrivendo su tutte le bande
gdal.Rasterize(
    ds_tif,
    shape,
    inverse=True,
    bands=list(range(1, ds_tif.RasterCount + 1)),
    burnValues=(0,) * ds_tif.RasterCount,
)
# Salvo il GeoTiff (alla fine di tutti i postprocessing)
gdal.Translate("output.png", ds_tif)
# NOTA: il file un memoria va chiuso manualmente. Si potrebbe fare una classe RAII
gdal.Unlink(tif_path)

edigiacomo commented 1 year ago

Il seguente codice permette di fare il clipping del PNG "georeferenziabile" generato da arkimaps usando uno shapefile qualsiasi caricato in memoria e senza appoggiarsi a file temporanei salvati su disco.

Dimenticavo: il ritaglio ha senso solo se il PNG è "georeferenziabile", quindi nel caso di tile e con flavour come emro_web (https://github.com/ARPA-SIMC/arkimaps/blob/v1.4-2/recipes/flavours/emro.yaml#L41). Negli altri casi non ha senso.

È facile identificare il flavour "tile" perché c'è il campo tile nella ricetta, ma questo non è vero nel caso di una singola immagine PNG georeferenziabile (come nel caso del flavour emro_web), perché il risultato è dato da una particolare configurazione dei parametri Magics. Si potrebbe far diventare questo flavour un "first class citizen" come i tile, estraendo la parte opportuna dalla ricetta e cablandola nel codice di arkimaps, analogamente ai tile.

spanezz commented 1 year ago

Intanto ho implementato una feature non richiesta come proof of concept di postprocessatori:

class Watermark(Postprocessor):
    """
    Write a string on the image.

    Parameters:

    * ``message``: text string to write
    * ``font``: name of a .ttf file to use as a font. The .ttf file needs to be
      found inside the static data directory
    * ``x``: horizontal coordinates (in pixel) of the beginning of the text. A
      negative value is the number of pixels from the right margin of the image
    * ``y``: vertical coordinates (in pixel) of the beginning of the text. A
      negative value is the number of pixels from the bottom margin of the image
    * ``anchor``: The text anchor alignment. See
      https://pillow.readthedocs.io/en/stable/handbook/text-anchors.html#text-anchors
    * ``size``: font size in pixels (default: 10)
    * ``color``: color name as defined in
      https://pillow.readthedocs.io/en/stable/reference/ImageColor.html#color-names
      Default: "#fff0"
    """

spanezz commented 1 year ago

Ho fatto push della prima implementazione del postprocessing cutshape (branch issue74):

 - name: emro_web
   postprocess:
     - type: cutshape
       shapefile: "er.geojson"

Non è un gran ché ottimizzata, però sembra funzionare.

Alcune note:

non ho fatto push di er.geojson nel repository perché è enorme. Può essere forse piú efficiente convertirlo in shapefile e aggiungere quello?
funziona solo quando in add_basemap si usano proiezioni epsg, altrimenti logga un warning e non fa nulla

E la nota piú importante di tutte: complimenti per il codice del prototipo di ritaglio :heart: Ho notato e apprezzato tante bellissime accortezze per renderlo, generico, succinto, e facile da integrare. È stato un piacere lavorarci!

edigiacomo commented 1 year ago

Grazie @spanezz, è un piacere sapere che è stato un piacere :)

non ho fatto push di er.geojson nel repository perché è enorme. Può essere forse piú efficiente convertirlo in shapefile e aggiungere quello?

Ho convertito il file da GeoJSON a Shapefile, lo trovi qui: er.zip. La dimensione passa da ~1.2MB a ~460K

funziona solo quando in add_basemap si usano proiezioni epsg, altrimenti logga un warning e non fa nulla

Per me va benissimo, se poi dobbiamo supportare altre proiezioni apriremo issue.

Grazie ancora!

spanezz commented 1 year ago

Perfetto, non solo prende meno spazio ma fa anche molto prima a caricarlo.

Ho fatto push del caricamento degli shape solo una volta per render script.

In recipes/flavours/emro.yaml ci sono postprocessing di esempio definiti: va bene finché il branch è un prototipo, ma vanno probabilmente tolti/sistemati con quello che effettivamente serve prima di un merge in main.

In particolare il flavour emro al momento non funziona perché richiede di mettere manualmente LiberationSans-Regular.ttf in recipes/static (non volevo aggiungere un ttf alla storia del repository solo per fare delle prove)

spanezz commented 1 year ago

Passo a voi per verifica

edigiacomo commented 1 year ago

Ho fatto delle prove con emro_web (commentando il postprocessing di emro) e direi che il ritaglio funziona perfettamente!

Per me è ok, puoi rimuovere il postprocess da emro e lasciare quello di emro_web.

edigiacomo commented 1 year ago

Ho sistemato i postprocessatori nei due flavour e il nome dello shapefile. Considero la issue risolta e faccio merge nel master.

ARPA-SIMC / arkimaps