Estrazione cartografia IGM non corretta tramite comandi GDAL

Saverio · February 1, 2025, 2:10am

L’obiettivo è quello di estrarre una parte di mappa IGM25000 dal geoportale cartografico italiano utilizzando GDAL_TRANSLATE

step 1
mi assicuro che i servizio funzioni e che le coordinate siano corrette tramite la seguente URL:

ottengo la visualizzazione della mappa desiderata in forma corretta come da seguente immagine:

step 2
tramite gdal_translate creazione file .xml che eviterà ogni volta di dare comandi con troppi parametri, memorizzando quelli costanti:

dalle capabilities del servizio, sotto la serie degli EPSG supportati, per l’ EPSG:4326 si ha:

Quindi BBOX= minx=35.5 miny=6.6 maxx=47.2 maxy=18.7 =35.5,6.6,47.2,18.7

gdal_translate -of WMS WMS:“http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&service=WMS&version=1.3.0&request=GetMap&layers=CB.IGM25000&format=image/jpeg&crs=EPSG:4326&bbox=35.5,6.6,47.2,18.7” “C:\D\prova\IGM25000-EPSG-4326.xml”

permette di creare il file IGM25000-EPSG-4326.xml che conterrà:

<GDAL_WMS>

1.3.0
MapServer Message
CB.IGM25000
EPSG:4326
image/jpeg
FALSE
xyXY

35.5
18.7
47.2
6.6
1038246226
1073741824

3
1024
1024
20
</GDAL_WMS>

step 3
estrazione mappa IGM tramite comando gdal_translate e utilizzo del file .xml creato sopra:

gdal_translate -projwin 45.8366237340226 9.47651291797857 45.8455710549957 9.46355987117446 -outsize 2000 2000 “C:\D\prova\IGM25000-EPSG-4326.xml” “C:\D\prova\IGM25000-geoportale-4326.tif”

Il file IGM25000-geoportale-4326.tif contiene la mappa IGM desiderata ma con tutte le tiles disordinate, ovvero senza una ricostruzione del mosaico corretto, che rende la mappa non leggibile.

step 4
fatte decine di prove invertendo l’ordine della sequenza delle latitudini e longitudini inserite sia sul comando gdal_translate allo step 3, sia all’interno del file .xml creato allo step 2.
riesco ad ottenere una rotazione della mappa estratta nel file .tif ma sempre con le tiles non riportate correttamente.
Evidenzio che seguendo la procedura sopra esposta ma utilizzando l’ EPSG:32632 con coordinate riferite al corrispondente SRS, ottengo una mappa estratta .tif CORRETTA!

Qualcuno può aiutarmi?
grazie

Saverio · February 1, 2025, 4:23pm

la URL allo step 1 è:
“wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&service=WMS&version=1.3.0&request=GetMap&layers=CB.IGM25000&format=image/jpeg&crs=EPSG:4326&bbox=45.8366237340226,9.46355987117446,45.8455710549957,9.47651291797857&width=2000&height=2000”

agiudiceandrea · February 1, 2025, 11:54pm

Ciao Saverio,
questa mailing list è relativa a QGIS, quindi penso che siamo un po’ fuori tema. Tuttavia penso anche che le tematiche siano attinenti.

Sarebbe utile se tu indicassi la versione di GDAL che stai utilizzando e come l’hai installata (per esempio, tramite OSGeo4W oppure diversamente).

A mio parere il problema è dovuto all’ordine in cui vengono espresse le coordinate.

Il comando da te indicato:
gdal_translate -of WMS WMS:“http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&service=WMS&version=1.3.0&request=GetMap&layers=CB.IGM25000&format=image/jpeg&crs=EPSG:4326&bbox=35.5,6.6,47.2,18.7” “C:\D\prova\IGM25000-EPSG-4326.xml”
nel mio sistema (Windows 10 - GDAL/OGR 3.10.1 - OSGeo4W) genera un file XML che contiene informazioni un po’ diverse e soprattutto contiene il tag <BBoxOrder>yxYX</BBoxOrder> che indica un ordine delle coordinate che è al contrario di quanto sembrerebbe essere presente nel tuo file XML (dove è indicato come xyXY, anche se testo non viene visualizzato correttamente perché non l’hai inglobato in triplici backtick):

<GDAL_WMS>
  <Service name="WMS">
    <Version>1.3.0</Version>
    <ServerUrl>http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&amp;service=WMS</ServerUrl>
    <Layers>CB.IGM25000</Layers>
    <CRS>EPSG:4326</CRS>
    <ImageFormat>image/jpeg</ImageFormat>
    <Transparent>FALSE</Transparent>
    <BBoxOrder>yxYX</BBoxOrder>
  </Service>
  <DataWindow>
    <UpperLeftX>6.6</UpperLeftX>
    <UpperLeftY>47.2</UpperLeftY>
    <LowerRightX>18.7</LowerRightX>
    <LowerRightY>35.5</LowerRightY>
    <SizeX>1073741824</SizeX>
    <SizeY>1038246226</SizeY>
  </DataWindow>
  <BandsCount>3</BandsCount>
  <BlockSizeX>1024</BlockSizeX>
  <BlockSizeY>1024</BlockSizeY>
  <OverviewCount>20</OverviewCount>
</GDAL_WMS>

Infatti l’EPSG:4326 ha l’ordine delle coordinate definito proprio come y,x (e non x,y come molti altri CRS).
Peraltro, nell’URL del WMS, tu stesso hai indicato le coordinate del bbox proprio come y,x,Y,X

Inoltre nel secondo comando che hai indicato:
gdal_translate -projwin 45.8366237340226 9.47651291797857 45.8455710549957 9.46355987117446 -outsize 2000 2000 “C:\D\prova\IGM25000-EPSG-4326.xml” “C:\D\prova\IGM25000-geoportale-4326.tif”
hai inserito in maniera errata le coordinate del parametro -projwin: per tale parametro, secondo il manuale di gdal_translate, le coordinate devono essere inserite nell’ordine <ulx> <uly> <lrx> <lry>, cioè upper left x, upper left y, lower right x, lower right y.

Nel mio sistema, dato il file XML generato precedentemente, il seguente comando (che è simile al quello indicato da te, tranne per il fatto che ha le coordinate del parametro -projwin nell’ordine corretto) genera un file GeoTIFF con i tiles ordinati correttamente:
gdal_translate -projwin 9.46355987117446 45.8455710549957 9.47651291797857 45.8366237340226 -outsize 2000 2000 C:\OSGeo4W_V2_B\IGM25000-EPSG-4326.xml C:\OSGeo4W_V2_B\IGM25000-geoportale-4326.tif

agiudiceandrea · February 2, 2025, 9:16am

Ciao Saverio,
mi sembra molto strano che gdal_translate produca due file XML diversi eseguendo lo stesso comando con la stessa identica versione della libreria GDAL/OGR.

Sarebbe utile se tu riprovassi ad eseguire lo stesso comando che hai indicato nel primo messaggio e verificassi il contenuto del file XML.

Comunque sia, ti segnalo che tale file XML, oltre a contenere il tag <BBoxOrder> con il valore errato xyXY (dovrebbe essere yxYX) contiene anche i tag <UpperLeftX>, <UpperLeftY>, <LowerRightX> e <LowerRightY> con valori errati (sono invertiti i valori di longitudine e latitudine) e conseguentemente anche i tag <BlockSizeX> e <BlockSizeY>.

Secondo me dovresti verificare se effettivamente tale problema occorre nel tuo sistema ed eventualmente risolverlo.

A presto.

Andrea

On another topic

Andrea infinite grazie per avermi dedicato del tempo e per aver risolto il problema.
L’errore nasce proprio dalla creazione errata del file .xml da parte del comando gdal_translate.
Anche io uso la versione 3.10.1 installata da OSGeo4W ma avevo fatto le prove anche con la 2.2.4 e con la 3.9.3 ma il formato del file .xml non cambiava rispetto a quello da me pubblicato e diverso dal tuo sia per il BBOX order che per la sequenza successiva delle coordinate per i 2 angoli del box da estrarre:
  <Service name="WMS">
    <Version>1.3.0</Version>
    <ServerUrl>http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&amp;service=WMS</ServerUrl>
    <Layers>CB.IGM25000</Layers>
    <CRS>EPSG:4326</CRS>
    <ImageFormat>image/jpeg</ImageFormat>
    <Transparent>FALSE</Transparent>
    <BBoxOrder>xyXY</BBoxOrder>
  </Service>
  <DataWindow>
    <UpperLeftX>35.5</UpperLeftX>
    <UpperLeftY>18.7</UpperLeftY>
    <LowerRightX>47.2</LowerRightX>
    <LowerRightY>6.6</LowerRightY>
    <SizeX>1038246226</SizeX>
    <SizeY>1073741824</SizeY>
  </DataWindow>
  <BandsCount>3</BandsCount>
  <BlockSizeX>1024</BlockSizeX>
  <BlockSizeY>1024</BlockSizeY>
  <OverviewCount>20</OverviewCount>
</GDAL_WMS>
Copiando la tua versione del file .xml funziona correttamente anche a me.
Io avevo dovuto inserire le coordinate al comando successivo nella sequenza non ordinata perché se la inserivo ordinata ottenevo un errore da parte di GDAL e non veniva prodotto l’output (perché non allineata al contenuto del file .xml).
Adesso che ho compreso il problema, anche se resta la diversità di comportamento della libreria gdal a parità di versione, cambierò manualmente il file .xml prodotto e per me il problema è risolto.
Di nuovo grazie e complimenti per la tua competenza!
Saverio

Saverio · February 2, 2025, 9:18am

Andrea infinite grazie per avermi dedicato del tempo e per aver risolto il problema.
L’errore nasce proprio dalla creazione errata del file .xml da parte del comando gdal_translate.
Anche io uso la versione 3.10.1 scaricata da trac.osgeo.org/osgeo4w/ ma avevo fatto le prove anche con la 2.2.4 e con la 3.9.3 ma il formato del file .xml non cambiava rispetto a quello da me pubblicato e diverso dal tuo sia per il BBOX order che per la sequenza successiva delle coordinate per i 2 angoli del box da estrarre:

  <Service name="WMS">
    <Version>1.3.0</Version>
    <ServerUrl>http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&amp;service=WMS</ServerUrl>
    <Layers>CB.IGM25000</Layers>
    <CRS>EPSG:4326</CRS>
    <ImageFormat>image/jpeg</ImageFormat>
    <Transparent>FALSE</Transparent>
    <BBoxOrder>xyXY</BBoxOrder>
  </Service>
  <DataWindow>
    <UpperLeftX>35.5</UpperLeftX>
    <UpperLeftY>18.7</UpperLeftY>
    <LowerRightX>47.2</LowerRightX>
    <LowerRightY>6.6</LowerRightY>
    <SizeX>1038246226</SizeX>
    <SizeY>1073741824</SizeY>
  </DataWindow>
  <BandsCount>3</BandsCount>
  <BlockSizeX>1024</BlockSizeX>
  <BlockSizeY>1024</BlockSizeY>
  <OverviewCount>20</OverviewCount>
</GDAL_WMS>

Copiando la tua versione del file .xml funziona correttamente anche a me.
Io avevo dovuto inserire le coordinate al comando successivo nella sequenza non ordinata perché se la inserivo ordinata ottenevo un errore da parte di GDAL e non veniva prodotto l’output (perché non allineata al contenuto del file .xml).
Adesso che ho compreso il problema, anche se resta la diversità di comportamento della libreria gdal a parità di versione, cambierò manualmente il file .xml prodotto e per me il problema è risolto. Poi se c’è qualche suggerimento per chiarire se esista qualche parametro per forzare Gdal nel formato corretto delle coordinate nella creazione del file .xml ben venga. O forse qualche parametro nel file di configurazione?
Potresti indicarmi cortesemente quale sia il blog più indicato per pubblicare post inerenti i comandi Gdal?
Di nuovo grazie e complimenti per la tua competenza! (io dopo diverse decine di prove non riuscivo ad uscirne fuori)
Saverio

agiudiceandrea · February 2, 2025, 5:50pm

Per verificare a cosa sia dovuto, dovresti indicare in dettaglio come esattamente hai installato GDAL (se hai usato l’installer di OSGeo4W, se hai usato la modalità Express o Advanced, quali pacchetti hai selezionato, …) e come esegui i comandi GDAL (per esempio, nel prompt dei comandi di Windows o nella Shell di OSGeo4W).

Potrebbe essere interessante anche eseguire il comando che genera il file xml aggiungendo l’opzione --debug ON in modo da verificare i messaggi di debug.

Saverio · February 3, 2025, 12:37am

Grazie Andrea per la tua ulteriore disponibilità.
Ho scaricato dal sito trac.osgeo.org/osgeo4w/ il software “OSGeo4W network installer”. Nella sua esecuzione sul mio PC ho utilizzato la modalità “advanced”, ho prima proceduto con l’opzione “uninstall” su tutti i moduli, compresi gli obsoleti, per fare pulizia delle precedenti versioni. Dopodiché ho anche cancellato la cartella dove era presente la precedente installazione,
Ho rilanciato nuovamente il software “OSGeo4W network installer” e questa volta ho selezionato la modalità “install” ma per sicurezza sull’ultima voce inerenti i moduli obsoleti ho di nuovo forzato “uninstall”.
L’installazione finisce correttamente con la frase “installazione di OSGEO4W completata con successo” anche se nel penultimo pannello era apparso il messaggio “package Shell - Shell.bat exit code 1” con la precisazione “questo non significa che i pacchetti interessati non funzioneranno correttamente”.
Nell’utilizzo dei comandi GDAL fino ad oggi non ho avuto problemi se non quest’ultimo sulla creazione del file .xml.
Di seguito l’output del comando gdal_translate con l’opzione --debug on per la creazione del file .xml.
Evidenzio che nonostante le operazioni di pulizia delle precedenti versioni come ho riportato sopra, nell’esecuzione dei comandi gdal continua ad uscire fuori il messaggio "ERROR 1: PROJ: proj_create_from_database: C:\OSGeo4W\share\proj\proj.db contains DATABASE.LAYOUT.VERSION.MINOR = 2 whereas a number >= 4 is expected. It comes from another PROJ installation. ma almeno apparentemente sembra un alert perché è sempre uscito con tutte le installazioni fatte per le varie versioni dalla 3 in poi ma i comandi gdal hanno sempre funzionato.

C:\OSGeo4W>gdalinfo --version
GDAL 3.10.1, released 2025/01/08

C:\OSGeo4W>
C:\OSGeo4W>gdal_translate -of WMS WMS:"http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&service=WMS&version=1.3.0&request=GetMap&layers=CB.IGM25000&format=image/jpeg&crs=EPSG:4326&bbox=35.5,6.6,47.2,18.7" "C:\D\prova\IGM25000-EPSG-4326.xml" --debug on
ERROR 1: PROJ: proj_create_from_database: C:\OSGeo4W\share\proj\proj.db contains DATABASE.LAYOUT.VERSION.MINOR = 2 whereas a number >= 4 is expected. It comes from another PROJ installation.
WMS: Opening WMS :
<GDAL_WMS>
  <Service name="WMS">
    <Version>1.3.0</Version>
    <ServerUrl>http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&amp;service=WMS</ServerUrl>
    <Layers>CB.IGM25000</Layers>
    <CRS>EPSG:4326</CRS>
    <ImageFormat>image/jpeg</ImageFormat>
    <Transparent>FALSE</Transparent>
    <BBoxOrder>xyXY</BBoxOrder>
  </Service>
  <DataWindow>
    <UpperLeftX>35.5</UpperLeftX>
    <UpperLeftY>18.7</UpperLeftY>
    <LowerRightX>47.2</LowerRightX>
    <LowerRightY>6.6</LowerRightY>
    <SizeX>1038246226</SizeX>
    <SizeY>1073741824</SizeY>
  </DataWindow>
  <BandsCount>3</BandsCount>
  <BlockSizeX>1024</BlockSizeX>
  <BlockSizeY>1024</BlockSizeY>
  <OverviewCount>20</OverviewCount>
</GDAL_WMS>

ERROR 1: PROJ: proj_create_from_database: C:\OSGeo4W\share\proj\proj.db contains DATABASE.LAYOUT.VERSION.MINOR = 2 whereas a number >= 4 is expected. It comes from another PROJ installation.
GDAL: GDALOpen(WMS:http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&service=WMS&version=1.3.0&request=GetMap&layers=CB.IGM25000&format=image/jpeg&crs=EPSG:4326&bbox=35.5,6.6,47.2,18.7, this=000001A522C49400) succeeds as WMS.
Input file size is 1038246226, 1073741824
GDAL: GDALOpen(C:\D\prova\IGM25000-EPSG-4326.xml, this=000001A522C49760) succeeds as WMS.
GDAL: GDALClose(C:\D\prova\IGM25000-EPSG-4326.xml, this=000001A522C49760)
WMS: Opening WMS :
<GDAL_WMS>
  <Service name="WMS">
    <Version>1.3.0</Version>
    <ServerUrl>http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&amp;service=WMS</ServerUrl>
    <Layers>CB.IGM25000</Layers>
    <CRS>EPSG:4326</CRS>
    <ImageFormat>image/jpeg</ImageFormat>
    <Transparent>FALSE</Transparent>
    <BBoxOrder>xyXY</BBoxOrder>
  </Service>
  <DataWindow>
    <UpperLeftX>35.5</UpperLeftX>
    <UpperLeftY>18.7</UpperLeftY>
    <LowerRightX>47.2</LowerRightX>
    <LowerRightY>6.6</LowerRightY>
    <SizeX>1038246226</SizeX>
    <SizeY>1073741824</SizeY>
  </DataWindow>
  <BandsCount>3</BandsCount>
  <BlockSizeX>1024</BlockSizeX>
  <BlockSizeY>1024</BlockSizeY>
  <OverviewCount>20</OverviewCount>
</GDAL_WMS>

GDAL: GDALOpen(WMS:http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&service=WMS&version=1.3.0&request=GetMap&layers=CB.IGM25000&format=image/jpeg&crs=EPSG:4326&bbox=35.5,6.6,47.2,18.7, this=000001A522C47FC0) succeeds as WMS.
GDAL: GDALClose(WMS:http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&service=WMS&version=1.3.0&request=GetMap&layers=CB.IGM25000&format=image/jpeg&crs=EPSG:4326&bbox=35.5,6.6,47.2,18.7, this=000001A522C47FC0)
GDAL: QuietDelete(C:\D\prova\IGM25000-EPSG-4326.xml) invoking Delete()
ERROR 1: PROJ: proj_create_from_database: C:\OSGeo4W\share\proj\proj.db contains DATABASE.LAYOUT.VERSION.MINOR = 2 whereas a number >= 4 is expected. It comes from another PROJ installation.
GDAL: GDALOpen(C:\D\prova\IGM25000-EPSG-4326.xml, this=000001A522C48320) succeeds as WMS.
GDAL: GDALClose(C:\D\prova\IGM25000-EPSG-4326.xml, this=000001A522C48320)
ERROR 1: PROJ: proj_create_from_database: C:\OSGeo4W\share\proj\proj.db contains DATABASE.LAYOUT.VERSION.MINOR = 2 whereas a number >= 4 is expected. It comes from another PROJ installation.
GDAL: GDALClose(C:\D\prova\IGM25000-EPSG-4326.xml, this=000001A522C490A0)
GDAL: GDALClose(WMS:http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_25000.map&service=WMS&version=1.3.0&request=GetMap&layers=CB.IGM25000&format=image/jpeg&crs=EPSG:4326&bbox=35.5,6.6,47.2,18.7, this=000001A522C49400)
GDAL: In GDALDestroy - unloading GDAL shared library.

C:\OSGeo4W> ```

I comandi GDAL li eseguo sempre nella shell di OSGEO4W lanciata con "osgeo4w.bat"

agiudiceandrea · February 3, 2025, 8:33am

Ciao Saverio,
il fatto che l’installer OSGeo4W ti abbia proposto di disinstallare i moduli obsoleti e il fatto che sia stato generato l’errore Shell.bat exit code 1 indicherebbe che la cartella in cui hai indicato di installare i pacchetti (presumibilmente C:\OSGeo4W\) conteneva dei file di una precedente installazione di vari anni fa e quindi non era stata completamente cancellata.

Inoltre il seguente messaggio di errore

indica che in C:\OSGeo4W\share\ è presente un file proj.db che, a giudicare dalla versione, è di almeno un anno fa.

La libreria GDAL usa la libreria PROJ per poter avere informazioni sui sistemi di riferimento e poterli gestire. A sua volta la libreria PROJ usa il file proj.db per ottenere informazioni sui sistemi di riferimento (tra cui, per esempio, il corretto ordine delle coordinate).

Quindi ti consiglio di rimuovere effettivamente e completamente la cartella C:\OSGeo4W\ prima di effettuare una installazione di GDAL oppure di indicare all’installer di OSGeo4W una diversa cartella (completamente vuota o non esistente) per l’installazione di GDAL.

A presto.

Andrea

Saverio · February 3, 2025, 4:02pm

Ciao Andrea, grazie per la risposta, come da mio precedente messaggio ho provveduto alla cancellazione delle precedenti versioni sia tramite l’installer di OSGEO4W (con l’opzione “uninstall” prima di procedere con l’opzione “install”), sia con la cancellazione manuale della directory C:\OSGEO4W.
Oggi ho provato anche con l’installazione su una directory diversa ma il risultato non cambia. Ho visto che in Internet ci sono diverse segnalazioni sull’ERROR 1 proj, ma nessuno indica una soluzione …
Saverio

agiudiceandrea · February 3, 2025, 5:02pm

Bisognerebbe leggere i log di installazione di OSGeo4W e l’elenco dei pacchetti installati e delle relative, ma questo sicuramente esula sia dal focus di questa mailing list sia dalla stessa possibilità che si possa risolvere il problema in tempi accettabili tramite scambio di messaggi in una mailing list, visto che la situazione da te descritta risulta incompatibile con la normale modalità in cui l’installazione dovrebbe procedere.

Mi rimane solo di consigliarti di creare un ticket su OSGeo4W e indicare dettagliatamente tutti i passaggi effettuati per l’installazione e allegando i log di installazione (in \var\log\) e il database dei pacchetti installati (\etc\setup\installed.db) e spiegando qual è il problema riscontrato (principalmente dovuto all’errore ERROR 1: PROJ: proj_create_from_database: C:\OSGeo4W\share\proj\proj.db contains DATABASE.LAYOUT.VERSION.MINOR = 2 whereas a number >= 4 is expected. It comes from another PROJ installation.)

Saverio · February 3, 2025, 5:51pm

Ok Andrea, grazie ancora per il tempo che mi hai dedicato e per la soluzione del problema principale che mi ha permesso comunque di poter estrarre oggi e in futuro la porzione di mappa IGM desiderata.
Saverio