[GIS] WMS e GetFeatureInfo

[pigreco]

Come posso scrapare/draping i dati di un WMS (http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/DTM_20M.map) avendo a disposizione un elenco di punti??

point.txt

[aborruso]

Caro @pigreco,
per queste cose puoi usare uno dei tesori che fornisce GDAL/OGR: gdallocationinfo.

Ti faccio un esempio su un WMS su cui ho lavorato da poco, che potrai adattare al caso:

gdallocationinfo "WMS:http://geoportal.asig.gov.al/service/wms?SERVICE=WMS&REQUEST=GetMap&SERVICE=WMS&VERSION=1.1.1&LAYERS=dem_2007:GridCoverage&FORMAT=image/png&SRS=EPSG:6870&BBOX=436435.5495804515,4389383.144576245,591627.5495804516,4725179.944576245" \
-xml -geoloc 499000 4550000  -b 1

Nell'URL è specificato il layer, il sistema di coordinate e il bounding box (l'estensione). Poi c'è il parametro per avere l'output in xml, poi con -geoloc 499000 4550000 dichiaro le coordinate del punto di cui voglio le info e con -b 1 indico la banda della raster di cui voglio le info (se rimuovi -b 1, puoi vedere tutte le bande e che info contengono).

In output

<Report pixel="200055709" line="560154335">
  <BandReport band="1">
    <LocationInfo>
      <wfs:FeatureCollection xmlns="http://www.opengis.net/wfs" xmlns:wfs="http://www.opengis.net/wfs" xmlns:dem_2007="/dem_2007" xmlns:gml="http://www.opengis.net/gml" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://geoportal.asig.gov.al/service/schemas/wfs/1.0.0/WFS-basic.xsd">
        <gml:boundedBy>
          <gml:null>unknown</gml:null>
        </gml:boundedBy>
        <gml:featureMember>
          <dem_2007:GridCoverage fid="">
            <dem_2007:ELEVATION>225.0</dem_2007:ELEVATION>
          </dem_2007:GridCoverage>
        </gml:featureMember>
      </wfs:FeatureCollection>
    </LocationInfo>
    <Value>255</Value>
  </BandReport>
</Report>

Qui il dato di interesse è la quota, che è in <dem_2007:ELEVATION>225.0</dem_2007:ELEVATION>.

Per estrarlo puoi usare scrape e fare una query XPATH. Siccome nella proprietà da estrarre ci sono i due punti ("dem_2007:ELEVATION") e i due punti sono speciali in un XML, li rimuoverò prima di fare la query.

gdallocationinfo "WMS:http://geoportal.asig.gov.al/service/wms?SERVICE=WMS&REQUEST=GetMap&SERVICE=WMS&VERSION=1.1.1&LAYERS=dem_2007:GridCoverage&FORMAT=image/png&SRS=EPSG:6870&BBOX=436435.5495804515,4389383.144576245,591627.5495804516,4725179.944576245" \
-xml -geoloc 499000 4550000  -b 1 | \
scrape -e '//elevation/text()'

Che ti darà 225.0.

Per farlo da un lista, non ti resta che creare un ciclo for loop in bash.

[aborruso]

@dataibi, @gpirrotta una query XPATH per estrarre dem_2007:ELEVATION (vedi XML presente qui) non può essere //dem_2007:ELEVATION, perché ci sono i due punti, le faccende di namespace, ecc..

Credo che il modo sia //*[name()='dem_2007:ELEVATION'], ma non credo che funzioni con oggetti basati su lxml.

C'è un altro modo, via XPATH?

Grazie

EDIT: ho trovato un modo

Arrivo al risultato con //*[contains(local-name(), "elevation")], ma non mi fa impazzire

[pigreco]

@aborruso
grazie per la dettagliata risposta e per il tempo dedicatomi. 👍

[aborruso]

@pigreco per me possiamo chiudere. Mi tengo la query XPATH bruttina