[Gfoss] cml2shp.py

Ciao a tutti,

Ho finalmente messo mano su un file CXF e uno CML, tutte e due
disponibile ugualmente dal portale per comuni del catasto. Ho capito
dal file (e non l'avevo capito dalla spec) che l'attributo e' gia'
logicamente associato alla geometria. Cosi il parsing e' molto piu'
facile che pensavo.

Allora mi ho messo a scrivere un piccolo parser per CML che e' xml e
cosi piu' facile che quello di CXF.

La versione attuale prende solo Bordi (incl. particelle e fabbricati).
Fin ora lo mette in una rappresentazione in memoria. Il prossimo passo
sara' di scrivere su un shape file usando la shapelib di python. Spero
di mandare una versione che include questo piu' tardi oggi.

Commenti molto benvenuti

saluti
-b

#######################################################################
## =======
## cml2shp
## =======
##
## Converts a subset of content from a file in
## Cadastral Markup Language (CML) [1] to a
## shape file
##
## [1] CML spec:
## http://www.agenziaterritorio.it/servizi/comunieistituzioni/
## ..fornitura_dati_catastali/specifica%20tecnica%20CML.doc
##
## Copyright: Comune di Grosseto
## Author: Bud P. Bruegger <bud@comune.grosseto.it>
## License: GPL (any version)
#######################################################################

import operator

import shapelib
import elementtree.ElementTree as ET

version = "0.1 14/12/2007"

#fileName = "testdata/E202_000100.CMF"
fileName = "testdata/vertisola.CMF"

class Feature(object):
    "a simple in memory rep of a feature"
    def __init__(self, bordoType, mappaID, particellaID, valenza,
esterconf, geom): self.bordoType = bordoType
        self.mappaID = mappaID
        self.particellaID = particellaID
        self.valenza = valenza
        self.esterconf = esterconf
        self.geom = geom
        #geometry is represented by a list of rings,
        #where every ring is a list of x/y tuples
        #the first ring is the main one,
        #the others are islands

#-------------------------------------
def parseFabbriPart(codbo):
    "parses a codbo of a fabbricato or particella"
    bordoType = "PARTICELLA"
    particellaID = codbo
    if codbo[-1] == '+':
        bordoType = "FABBRICATO"
        particellaID = codbo[:-1]
    return bordoType, particellaID
        
def codboParse(codbo, mapName):
    "parses a codbo"
    #returns: bordoType, particellaID
    codbo = codbo.strip()
    mappaID = mapName
    if codbo == mapName:
        bordoType = "MAPPA"
        particellaID = ""
    elif codbo == "STRADA":
        bordoType = "STRADA"
        particellaID = ""
    elif codbo == "ACQUA":
        bordoType = "ACQUA"
        particellaID = ""
    else:
        bordoType, particellaID = parseFabbriPart(codbo)
    return bordoType, particellaID

def vertConsume(nVerts, vertList):
    return vertList[:nVerts], vertList[nVerts:]

def coordParse(nVert, vertIsolaList, coordStr):
    noIsolaVerts = reduce(operator.add, vertIsolaList, 0)
    vertexList = [subStr.split(',') for subStr in coordStr.split()]
    geom =
    mainRing, remainingVerts = vertConsume((nVert - noIsolaVerts),
vertexList) geom.append(mainRing)
    for noIslandVerts in vertIsolaList:
        islandRing, remainingVerts = vertConsume(noIslandVerts,
remainingVerts) geom.append(islandRing)
    return geom

#-- main parsing ---------------------
mapTree = ET.parse(fileName)
mapRoot = mapTree.getroot()
mapID = mapRoot.find("INFOMAPPA").get("nome")
featureList =

#-- process BORDO elements ------------------
for border in mapRoot.getiterator("BORDO"):
    bordoType, particellaID = codboParse(border.get("codbo"), mapID)
    valenza = border.get("valenza")
    esterconf = border.get("esterconf")
    gbordo = border.find("GBORDO")
    #nIsole = gbordo.get("n.isole")
    nVert = int(gbordo.get("n.vert"))
    vertIsolaList = map(int, [i.text for i in gbordo.getiterator
("VERTISOLA")]) coordStr = gbordo.find("COORD").text
    geom = coordParse(nVert, vertIsolaList, coordStr)
    featureList.append(Feature(bordoType, mapID, particellaID,
                               valenza, esterconf, geom))

--
Bud P. Bruegger, Ph.D. +39-0564-488577 (voice), -21139 (fax)
   European Chair, Global Collaboration Forum on eID
   Chair, Porvoo Subgroup on collab. govs/operating systems
   Leader of the Permanent eID Status Observatory (PESO) project
Servizio Elaborazione Dati e-mail: bud@comune.grosseto.it
Comune di Grosseto jabber: bud@jabber.no
Via Ginori, 43 http://www.comune.grosseto.it/
58100 Grosseto (Tuscany, Italy)
http://www.comune.grosseto.it/interopEID/

Ciao Bud,

complimenti davvero per il lavoro.
Potremo svilupparlo come plugin per Qgis in modo da avere lo shape già caricato nella mappa..

Se ti va posso darti una mano.

Luca

Bud P. Bruegger ha scritto:

Ciao a tutti,

Ho finalmente messo mano su un file CXF e uno CML, tutte e due
disponibile ugualmente dal portale per comuni del catasto. Ho capito
dal file (e non l'avevo capito dalla spec) che l'attributo e' gia'
logicamente associato alla geometria. Cosi il parsing e' molto piu'
facile che pensavo.

Allora mi ho messo a scrivere un piccolo parser per CML che e' xml e
cosi piu' facile che quello di CXF.

La versione attuale prende solo Bordi (incl. particelle e fabbricati).
Fin ora lo mette in una rappresentazione in memoria. Il prossimo passo
sara' di scrivere su un shape file usando la shapelib di python. Spero
di mandare una versione che include questo piu' tardi oggi.

Commenti molto benvenuti

saluti
-b

#######################################################################
## =======
## cml2shp
## =======
##
## Converts a subset of content from a file in
## Cadastral Markup Language (CML) [1] to a
## shape file
##
## [1] CML spec:
## http://www.agenziaterritorio.it/servizi/comunieistituzioni/
## ..fornitura_dati_catastali/specifica%20tecnica%20CML.doc
##
## Copyright: Comune di Grosseto
## Author: Bud P. Bruegger <bud@comune.grosseto.it>
## License: GPL (any version)
#######################################################################

import operator

import shapelib
import elementtree.ElementTree as ET

version = "0.1 14/12/2007"

#fileName = "testdata/E202_000100.CMF"
fileName = "testdata/vertisola.CMF"

class Feature(object):
    "a simple in memory rep of a feature"
    def __init__(self, bordoType, mappaID, particellaID, valenza,
esterconf, geom): self.bordoType = bordoType
        self.mappaID = mappaID
        self.particellaID = particellaID
        self.valenza = valenza
        self.esterconf = esterconf
        self.geom = geom
        #geometry is represented by a list of rings,
        #where every ring is a list of x/y tuples
        #the first ring is the main one,
        #the others are islands

#-------------------------------------
def parseFabbriPart(codbo):
    "parses a codbo of a fabbricato or particella"
    bordoType = "PARTICELLA"
    particellaID = codbo
    if codbo[-1] == '+':
        bordoType = "FABBRICATO"
        particellaID = codbo[:-1]
    return bordoType, particellaID
        def codboParse(codbo, mapName):
    "parses a codbo"
    #returns: bordoType, particellaID
    codbo = codbo.strip()
    mappaID = mapName
    if codbo == mapName:
        bordoType = "MAPPA"
        particellaID = ""
    elif codbo == "STRADA":
        bordoType = "STRADA"
        particellaID = ""
    elif codbo == "ACQUA":
        bordoType = "ACQUA"
        particellaID = ""
    else:
        bordoType, particellaID = parseFabbriPart(codbo)
    return bordoType, particellaID

def vertConsume(nVerts, vertList):
    return vertList[:nVerts], vertList[nVerts:]

def coordParse(nVert, vertIsolaList, coordStr):
    noIsolaVerts = reduce(operator.add, vertIsolaList, 0)
    vertexList = [subStr.split(',') for subStr in coordStr.split()]
    geom =
    mainRing, remainingVerts = vertConsume((nVert - noIsolaVerts),
vertexList) geom.append(mainRing)
    for noIslandVerts in vertIsolaList:
        islandRing, remainingVerts = vertConsume(noIslandVerts,
remainingVerts) geom.append(islandRing)
    return geom

#-- main parsing ---------------------
mapTree = ET.parse(fileName)
mapRoot = mapTree.getroot()
mapID = mapRoot.find("INFOMAPPA").get("nome")
featureList =

#-- process BORDO elements ------------------
for border in mapRoot.getiterator("BORDO"):
    bordoType, particellaID = codboParse(border.get("codbo"), mapID)
    valenza = border.get("valenza")
    esterconf = border.get("esterconf")
    gbordo = border.find("GBORDO")
    #nIsole = gbordo.get("n.isole")
    nVert = int(gbordo.get("n.vert"))
    vertIsolaList = map(int, [i.text for i in gbordo.getiterator
("VERTISOLA")]) coordStr = gbordo.find("COORD").text
    geom = coordParse(nVert, vertIsolaList, coordStr)
    featureList.append(Feature(bordoType, mapID, particellaID,
                               valenza, esterconf, geom))

Luca Casagrande ha scritto:

Potremo svilupparlo come plugin per Qgis in modo da avere lo shape già
caricato nella mappa..

Se ti va posso darti una mano.

Bravi!!
pc
--
Paolo Cavallini, see: http://www.faunalia.it/pc

Visto che siamo in tanti, potremo organizzare un gruppo di lavoro sull'argomento in
modo da dividerci i compiti e velocizzare lo sviluppo.

Metto su una pagina nel Wiki?

Ciao
Luca

Ciao

On Fri, 14 Dec 2007 13:13:10 +0100
Luca Casagrande <luca.casagrande@gmail.com> wrote:

Ciao Bud,

complimenti davvero per il lavoro.

Come hai visto era una cosa veloce..

Potremo svilupparlo come plugin per Qgis in modo da avere lo shape già
caricato nella mappa..

Non so ancora niente di Qgis ma suona molto bene.

Se ti va posso darti una mano.

Ho alcuni problemi ancora con il backend per scrivere shape files. Ho
solo guardato un esempio del uso di shapelib senza leggere la spec...

Che suggede e' che shpdump mi fa vedere tutto ok sembra a prima vista.
QGIS non mi fa vedere niente, e usa la scala di "degrees" assumendo
lat/long penso. Cosi una possibilita' e' che assumendo coordinate
geografiche non fa vedere niente oppure non ho fatto i "rings" nella
seguenza giusta (ho sparato senza vedere...).

C'e' qualcuno che ha idea che potrebbe essere spagliato?

Il codice fin ora e' allegato..

saluti
-b

cmf2shp.py (4.4 KB)

Purtroppo non conosco python e quindi mi e' difficile seguire il flusso del codice che hai scritto.

Comunque,

senz'altro devi rivedere i nomi dei campi DBF perche' la specifica DBF ammette al max 10 caratteri per il nome di un campo di attributi,
e di conseguenza anche per lo shapefile e' uguale.

e quando definisci il campo

dbf.add_field("particellaID", dbflib.FTString, 20, 0)

impegni 12 caratteri.

Per il resto mi pare di capire che usi la shapelib,
che pero' , innestata dentro la sintassi python mi e' completamente irriconoscibile.

E quindi non posso darti nessun feedback utile.

Posso solo suggerirti, se non lo fa' in automatico python, di distruggere l'oggetto SHPObject dopo averlo scritto sul file,
nella shapelib si usa "SHPDestroyObject()".
Altrimenti, se generi dei files con molti oggetti, andrai rapidamente ad esaurire le risorse della macchina.

Andrea.

On Fri, 14 Dec 2007 15:59:23 +0100, Bud P. Bruegger wrote:

Ciao

On Fri, 14 Dec 2007 13:13:10 +0100
Luca Casagrande <luca.casagrande@gmail.com> wrote:

Ciao Bud,

complimenti davvero per il lavoro.

Come hai visto era una cosa veloce..

Potremo svilupparlo come plugin per Qgis in modo da avere lo shape già
caricato nella mappa..

Non so ancora niente di Qgis ma suona molto bene.

Se ti va posso darti una mano.

Ho alcuni problemi ancora con il backend per scrivere shape files. Ho
solo guardato un esempio del uso di shapelib senza leggere la spec...

Che suggede e' che shpdump mi fa vedere tutto ok sembra a prima vista.
QGIS non mi fa vedere niente, e usa la scala di "degrees" assumendo
lat/long penso. Cosi una possibilita' e' che assumendo coordinate
geografiche non fa vedere niente oppure non ho fatto i "rings" nella
seguenza giusta (ho sparato senza vedere...).

C'e' qualcuno che ha idea che potrebbe essere spagliato?

Il codice fin ora e' allegato..

saluti
-b

Bud P. Bruegger ha scritto:

Ho alcuni problemi ancora con il backend per scrivere shape files. Ho
solo guardato un esempio del uso di shapelib senza leggere la spec...

Che suggede e' che shpdump mi fa vedere tutto ok sembra a prima vista.
QGIS non mi fa vedere niente, e usa la scala di "degrees" assumendo
lat/long penso. Cosi una possibilita' e' che assumendo coordinate
geografiche non fa vedere niente oppure non ho fatto i "rings" nella
seguenza giusta (ho sparato senza vedere...).

C'e' qualcuno che ha idea che potrebbe essere spagliato?

Il codice fin ora e' allegato..

Complimenti Bud! Mi sembra un ottimo inizio.
Per quanto riguarda il problema di Qgis, per avere una rappresentazione
più o meno congruente, ti conviene assegnare la proiezione epsg:26591
per la vista (anche se deprecato dall'epsg dal 2003, rappresenta pur
sempre il ns sistema cartografico nazionale... senza aggiungere il
parametro +towgs84 che definisce la trasformazione verso WGS84),
abilitare la proiezione al volo e creare infine una proiezione
personalizzata da assegnare ai catastali del tipo:
+proj=cass +lat_0=... +lon_0=... +x_0=... +y_0=... +ellps=intl +pm=rome
+units=m +no_defs
dove:
- +lat_0 e +lon_0 sono le coordinate geografiche del centro di
emanazione nel datum della vista (epsg:4806);
- +x_0 e +y_0 le falsi origini del tuo sistema catastale locale
(generalmente sono pari a 0 entrambi);
- +ellps=bessel poichè Bessel è l'ellissoide di Cassini-Soldner.
Questa procedura ti consente, con le dovute approssimazioni del caso, di
proiettare al volo i catastali in Gauss-Boaga Roma40 Fuso Ovest.
L'errore relativo seguendo tale procedura può essere anche di qualche
metro, in particolare aumenta all'aumentare della distanza dal meridiano
passante per il centro di emanazione.
Analogo è il caso del fuso Est (come nel mio caso), basta semplicemente
utilizzare 26592 per la vista ed il resto è lo stesso. Meglio di niente!
Per Luca: la pagina nel wiki è d'obbligo! :wink:

Ciao
Antonio

più o meno congruente, ti conviene assegnare la proiezione epsg:26591
per la vista (anche se deprecato dall'epsg dal 2003, rappresenta pur

EPSG ha deprecato (e tolto) i codici 26591 e 26592, dice di usare
3003 e 3004 al loro posto. Che io sappia non ci sono
controindicazioni.

Ho messo queste info anche su:
http://it.wikipedia.org/wiki/Proiezione_di_Gauss-Boaga

--
Niccolo Rigacci
Firenze - Italy

Salve a tutti.
Navigando un pò oggi mi sono imbattuto nell'implementazione del WPS in Java (che vorrei testare)
http://52north.org/index.php?option=com_projects&task=showProject&id=21&Itemid=127
e del relativo plugin di connessione per uDIG:
https://52north.org/twiki/bin/view/Processing/52nUdigWPSClient

Che si dice per gvSIG?

Ciao

Grazie Andrea,

certo un errore da correggere.

Ho anche riletto la spec di CML e ho notato che ci possono essere
moltipli elementi di "BORDO". Una altra cosa da correggere..

La parte di parsing dovrebbe funzionare bene, ma non e' testato. La
parte di backend e' molto rozzo--ho provato (ancora senza successo) di
far visualizzare qualcosa al volo :wink:

Devo ancora torvare un po di documentazione di shapelib in python e/o
esempi..

saluti e grazie
-b

On Sat, 15 Dec 2007 10:11:52 +0100
"Andrea P." <cerebrogis@ipergeo.org> wrote:

Purtroppo non conosco python e quindi mi e' difficile seguire il flusso del codice che hai scritto.

Comunque,

senz'altro devi rivedere i nomi dei campi DBF perche' la specifica DBF ammette al max 10 caratteri per il nome di un campo di attributi,
e di conseguenza anche per lo shapefile e' uguale.

e quando definisci il campo
>dbf.add_field("particellaID", dbflib.FTString, 20, 0)
impegni 12 caratteri.

Per il resto mi pare di capire che usi la shapelib,
che pero' , innestata dentro la sintassi python mi e' completamente irriconoscibile.

E quindi non posso darti nessun feedback utile.

Posso solo suggerirti, se non lo fa' in automatico python, di distruggere l'oggetto SHPObject dopo averlo scritto sul file,
nella shapelib si usa "SHPDestroyObject()".
Altrimenti, se generi dei files con molti oggetti, andrai rapidamente ad esaurire le risorse della macchina.

Andrea.

On Fri, 14 Dec 2007 15:59:23 +0100, Bud P. Bruegger wrote:

>Ciao

>On Fri, 14 Dec 2007 13:13:10 +0100
>Luca Casagrande <luca.casagrande@gmail.com> wrote:

>> Ciao Bud,
>>
>> complimenti davvero per il lavoro.

>Come hai visto era una cosa veloce..

>> Potremo svilupparlo come plugin per Qgis in modo da avere lo shape già
>> caricato nella mappa..

>Non so ancora niente di Qgis ma suona molto bene.

>> Se ti va posso darti una mano.

>Ho alcuni problemi ancora con il backend per scrivere shape files. Ho
>solo guardato un esempio del uso di shapelib senza leggere la spec...

>Che suggede e' che shpdump mi fa vedere tutto ok sembra a prima vista.
>QGIS non mi fa vedere niente, e usa la scala di "degrees" assumendo
>lat/long penso. Cosi una possibilita' e' che assumendo coordinate
>geografiche non fa vedere niente oppure non ho fatto i "rings" nella
>seguenza giusta (ho sparato senza vedere...).

>C'e' qualcuno che ha idea che potrebbe essere spagliato?

>Il codice fin ora e' allegato..

>saluti
>-b

_______________________________________________
Prenota la tua maglietta GFOSS.it:
http://wiki.gfoss.it/index.php/Gadgets
Iscriviti all'associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@faunalia.com
http://www.faunalia.com/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell'Associazione GFOSS.it.

--
Bud P. Bruegger, Ph.D. +39-0564-488577 (voice), -21139 (fax)
   European Chair, Global Collaboration Forum on eID
   Chair, Porvoo Subgroup on collab. govs/operating systems
   Leader of the Permanent eID Status Observatory (PESO) project
Servizio Elaborazione Dati e-mail: bud@comune.grosseto.it
Comune di Grosseto jabber: bud@jabber.no
Via Ginori, 43 http://www.comune.grosseto.it/
58100 Grosseto (Tuscany, Italy)
http://www.comune.grosseto.it/interopEID/

Ing. Fabio D'Ovidio ha scritto:

Salve a tutti.
Navigando un pò oggi mi sono imbattuto nell'implementazione del WPS in Java (che vorrei testare)
http://52north.org/index.php?option=com_projects&task=showProject&id=21&Itemid=127
e del relativo plugin di connessione per uDIG:
https://52north.org/twiki/bin/view/Processing/52nUdigWPSClient

Che si dice per gvSIG?

Fabio, in giro non mi è parso di vedere tracce di implementazioni WPS
per gvSIG, ma in ogni caso penso che il riuso del plugin di uDig in gvSIG
sia sicuramente possibile, visto che utilizzano molte librerie/driver in
comune, apportando le dovute modifiche del caso.

Ciao
Antonio

On Dec 19, 2007, at 10:52 AM, Antonio Falciano wrote:

Ing. Fabio D'Ovidio ha scritto:

Salve a tutti.
Navigando un pò oggi mi sono imbattuto nell'implementazione del WPS in
Java (che vorrei testare)
http://52north.org/index.php?option=com_projects&task=showProject&id=21&Itemid=127
e del relativo plugin di connessione per uDIG:
https://52north.org/twiki/bin/view/Processing/52nUdigWPSClient

Che si dice per gvSIG?

Fabio, in giro non mi è parso di vedere tracce di implementazioni WPS
per gvSIG, ma in ogni caso penso che il riuso del plugin di uDig in gvSIG
sia sicuramente possibile, visto che utilizzano molte librerie/driver in
comune, apportando le dovute modifiche del caso.

Hmmm... la vedo piu' facile partendo dallo stesso plugin che i 52north hanno fatto per jump, visto che credo abbiano piu' librerie in comune.

Ciao
Andrea

Ciao
Antonio

_______________________________________________
Prenota la tua maglietta GFOSS.it:
http://wiki.gfoss.it/index.php/Gadgets
Iscriviti all'associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@faunalia.com
http://www.faunalia.com/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell'Associazione GFOSS.it.

Andrea Antonello ha scritto:

On Dec 19, 2007, at 10:52 AM, Antonio Falciano wrote:

Ing. Fabio D'Ovidio ha scritto:

Salve a tutti.
Navigando un pò oggi mi sono imbattuto nell'implementazione del WPS in
Java (che vorrei testare)
http://52north.org/index.php?option=com_projects&task=showProject&id=21&Itemid=127
e del relativo plugin di connessione per uDIG:
https://52north.org/twiki/bin/view/Processing/52nUdigWPSClient

Che si dice per gvSIG?

Fabio, in giro non mi è parso di vedere tracce di implementazioni WPS
per gvSIG, ma in ogni caso penso che il riuso del plugin di uDig in gvSIG
sia sicuramente possibile, visto che utilizzano molte librerie/ driver in
comune, apportando le dovute modifiche del caso.

Hmmm... la vedo piu' facile partendo dallo stesso plugin che i 52north hanno fatto per jump, visto che credo abbiano piu' librerie in comune.

Meglio ancora! Massimo risultato con il minimo sforzo! :wink:

Ciao
Antonio