[Gfoss] verifica senso digitalizzazione geometria

Salve,
ho un problema con alcune geometrie areali, avrei la necessità di verificare il senso di digitalizzazione (orario:interno poligono pieno; antiorario: interno poligono vuoto);
esiste quest possibilità in Qgis?

Grazie mille,
marco


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi geografici
applicazioni GIS e webGIS

Ciao Marco,

non sono ancora praticissimo di QGIS, ma se usi la nuova simbologia,
puoi aggiungere allo stile poligono un nuovo layer di tipo "Line
decoration" e posizionarlo sopra il "Simple fill":

Layer properties | Style | New symbology | Change...

clicca su +, seleziona il layer aggiunto, e sulla combo Layer Tipe
selezioni Line decoration.

Correntemente tra i line decoration io posso selezionare solo una
freccina puntata verso destra, che immagino essere orientata lungo il
verso delle geometrie. Nel qual caso sei a posto. Puoi cortesemente
verificarlo?

Ne approfitto per aggiungere che questa possibilità di comporre la
simbologia in QGIS è molto bella anche se devo ancora capire bene come
funziona, in particolare come aggiungere nuovi "mattoncini" oltre che
comporre quelli preesistenti.

E' possibile ad esempio ripetere la freccina a intervalli regolari lungo
la linea? E avere la freccina girata in senso contrario alla geometria?
Mi sarebbe utilissima per tematizzare i sensi di marcia dei grafi
stradali...

Ciao a grazie a tutti

Sig

Il giorno mer, 15/06/2011 alle 11.00 +0200, marco zanieri ha scritto:

Salve,
ho un problema con alcune geometrie areali, avrei la necessità di
verificare il senso di digitalizzazione (orario:interno poligono
pieno; antiorario: interno poligono vuoto);
esiste quest possibilità in Qgis?

Grazie mille,
marco

--
            dott. Marco Zanieri
   e-mail: marcozanieri@gmail.com

           cartografia tematica
          banche dati territoriali
     sistemi informativi geografici
      applicazioni GIS e webGIS

_______________________________________________
Iscriviti all'associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell'Associazione GFOSS.it.
518 iscritti al 3.6.2011

Ciao Marco,

ho fatto ulteriori esperimenti, funziona meglio se aggiungi un layer di
tipo "Marker Line", e selezioni il simbolo ">".

Per rispondermi, si può applicare una rotazione di 180° per fare in modo
che la freccia punti al verso contrario della linea. E si può ripetere
il simbolo a intervalli regolari così come posizionarlo sui vertici o
nel punto medio.

Il problema è che lungo i lati dei poligoni adiacenti le frecce si
sovrappongono e non sono di lettura facile. Quindi in questo caso ti
consiglio il posizionamento sul punto medio.

Per evitare che il poligono adiacente si sovrapponga al tuo simbolo,
attiva i "Symbol Levels"

Ciao

Sig

Il giorno mer, 15/06/2011 alle 11.00 +0200, marco zanieri ha scritto:

Salve,
ho un problema con alcune geometrie areali, avrei la necessità di
verificare il senso di digitalizzazione (orario:interno poligono
pieno; antiorario: interno poligono vuoto);
esiste quest possibilità in Qgis?

Grazie mille,
marco

--
            dott. Marco Zanieri
   e-mail: marcozanieri@gmail.com

           cartografia tematica
          banche dati territoriali
     sistemi informativi geografici
      applicazioni GIS e webGIS

_______________________________________________
Iscriviti all'associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell'Associazione GFOSS.it.
518 iscritti al 3.6.2011

Dall’API di QGis non trovo un metodo per determinare l’orientamento della sequenza di coordinate… o se si tratta di un’isola. Un modo grezzo ma, credo, efficace è lanciare questo script nella console di Python dentro QGis, nel quale faccio un intersect tra la geometria e il suo centroide. Se fossse un’isola dovrebbe tornarmi False… o sbaglio?

iface = qgis.utils.iface
lyr = iface.activeLayer()
prov = lyr.dataProvider()
attrlist = prov.attributeIndexes()
prov.select(attrlist)
feat = QgsFeature()
for i in range(lyr.featureCount()):
prov.nextFeature(feat)
geom = feat.geometry()
cent = geom.centroid()
dentro = geom.intersects(cent)
if (dentro):
attributes = feat.attributeMap()
print ‘La feature %s contiene una geometria piena’ % attributes[0].toInt()[0]
else:
print ‘->>> La feature %s sembra contenere un’isola’ % attributes[0].toInt()[0]

Giovanni

Il giorno 15 giugno 2011 11:00, marco zanieri <marcozanieri@gmail.com> ha scritto:

Salve,
ho un problema con alcune geometrie areali, avrei la necessità di verificare il senso di digitalizzazione (orario:interno poligono pieno; antiorario: interno poligono vuoto);
esiste quest possibilità in Qgis?

Grazie mille,
marco


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi geografici
applicazioni GIS e webGIS


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011

Il test del centroide è un test del pirla, scusate, perché il centroide può cadere anche fuori dalla geometria…
Posterò una domanda nell ml di sviluppo ma intanto lo chiedo anche qua: come si può sapere l’orientamento di una geometria tramite le api di Qgis?

giovanni

Il giorno 15 giugno 2011 12:36, G. Allegri <giohappy@gmail.com> ha scritto:

Dall’API di QGis non trovo un metodo per determinare l’orientamento della sequenza di coordinate… o se si tratta di un’isola. Un modo grezzo ma, credo, efficace è lanciare questo script nella console di Python dentro QGis, nel quale faccio un intersect tra la geometria e il suo centroide. Se fossse un’isola dovrebbe tornarmi False… o sbaglio?

iface = qgis.utils.iface
lyr = iface.activeLayer()
prov = lyr.dataProvider()
attrlist = prov.attributeIndexes()
prov.select(attrlist)
feat = QgsFeature()
for i in range(lyr.featureCount()):
prov.nextFeature(feat)
geom = feat.geometry()
cent = geom.centroid()
dentro = geom.intersects(cent)
if (dentro):
attributes = feat.attributeMap()
print ‘La feature %s contiene una geometria piena’ % attributes[0].toInt()[0]
else:
print ‘->>> La feature %s sembra contenere un’isola’ % attributes[0].toInt()[0]

Giovanni

Il giorno 15 giugno 2011 11:00, marco zanieri <marcozanieri@gmail.com> ha scritto:

Salve,
ho un problema con alcune geometrie areali, avrei la necessità di verificare il senso di digitalizzazione (orario:interno poligono pieno; antiorario: interno poligono vuoto);
esiste quest possibilità in Qgis?

Grazie mille,
marco


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi geografici
applicazioni GIS e webGIS


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011

ho provato con la tematizzazione delle coordinate e si riesce a comprendere l’orientamento…grazie mille Luca; credo comunque che possa risultare molto utile la possibilità di individuare la sequenza delle coordinate di una geometria e comunque una procedura automatica che possa riconoscere la presenza di eventuali isole (pieni e vuoti)…
grazie a tutti

Il giorno 15 giugno 2011 12:49, G. Allegri <giohappy@gmail.com> ha scritto:

Il test del centroide è un test del pirla, scusate, perché il centroide può cadere anche fuori dalla geometria…
Posterò una domanda nell ml di sviluppo ma intanto lo chiedo anche qua: come si può sapere l’orientamento di una geometria tramite le api di Qgis?

giovanni

Il giorno 15 giugno 2011 12:36, G. Allegri <giohappy@gmail.com> ha scritto:

Dall’API di QGis non trovo un metodo per determinare l’orientamento della sequenza di coordinate… o se si tratta di un’isola. Un modo grezzo ma, credo, efficace è lanciare questo script nella console di Python dentro QGis, nel quale faccio un intersect tra la geometria e il suo centroide. Se fossse un’isola dovrebbe tornarmi False… o sbaglio?

iface = qgis.utils.iface
lyr = iface.activeLayer()
prov = lyr.dataProvider()
attrlist = prov.attributeIndexes()
prov.select(attrlist)
feat = QgsFeature()
for i in range(lyr.featureCount()):
prov.nextFeature(feat)
geom = feat.geometry()
cent = geom.centroid()
dentro = geom.intersects(cent)
if (dentro):
attributes = feat.attributeMap()
print ‘La feature %s contiene una geometria piena’ % attributes[0].toInt()[0]
else:
print ‘->>> La feature %s sembra contenere un’isola’ % attributes[0].toInt()[0]

Giovanni

Il giorno 15 giugno 2011 11:00, marco zanieri <marcozanieri@gmail.com> ha scritto:

Salve,
ho un problema con alcune geometrie areali, avrei la necessità di verificare il senso di digitalizzazione (orario:interno poligono pieno; antiorario: interno poligono vuoto);
esiste quest possibilità in Qgis?

Grazie mille,
marco


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi geografici
applicazioni GIS e webGIS


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi geografici
applicazioni GIS e webGIS

Ciao a tutti,

in teoria si fa così:

- estraggo il linearring del poligono;
- ne genero il centroide
- determino la proiezione e la segmentazione del centroide sulla
polilinea che corrisponde al linearring.

Se il centroide cade in "destra geometrica", allora l'arco è
digitalizzato in senso orario; altrimenti, in senso antiorario.

Per proiezione e segmentazione, intendo ad esempio le funzioni LRS di
PostGIS (ST_Line_Locate_Point() per la segmentazione, ma non trovo
quella che mi dà il lato). Ci sono funzioni analoghe in JTS, quindi in
GEOS e quindi in python. Secoli fa scrissi una classe Java per gestire
la segmentazione usando JTS, se ti serve te la spedisco, ma spero che
qualcuno abbia già implementato una funzione ad alto livello che dato un
punto e una polilinea mi dice a che distanza, da che lato e con che
proiezione questo cade su quella. Le funzioni JTS lavorano a livello di
punto e segmento (2 punti), quindi devi fare un ciclo per lavorare su
una linestring.

Fammi sapere se trovi qualcosa, mi interessa molto (soprattutto mi piace
l'idea di non dover fare il porting di quanto già implementato anni fa e
potermi così dedicare a cose nuove)

Sig

Il giorno mer, 15/06/2011 alle 13.00 +0200, marco zanieri ha scritto:

ho provato con la tematizzazione delle coordinate e si riesce a
comprendere l'orientamento...grazie mille Luca; credo comunque che
possa risultare molto utile la possibilità di individuare la sequenza
delle coordinate di una geometria e comunque una procedura automatica
che possa riconoscere la presenza di eventuali isole (pieni e
vuoti)...
grazie a tutti

Il giorno 15 giugno 2011 12:49, G. Allegri <giohappy@gmail.com> ha
scritto:
        Il test del centroide è un test del pirla, scusate, perché il
        centroide può cadere anche fuori dalla geometria...
        Posterò una domanda nell ml di sviluppo ma intanto lo chiedo
        anche qua: come si può sapere l'orientamento di una geometria
        tramite le api di Qgis?
        
        giovanni
        
        Il giorno 15 giugno 2011 12:36, G. Allegri
        <giohappy@gmail.com> ha scritto:
        
                Dall'API di QGis non trovo un metodo per determinare
                l'orientamento della sequenza di coordinate... o se si
                tratta di un'isola. Un modo grezzo ma, credo, efficace
                è lanciare questo script nella console di Python
                dentro QGis, nel quale faccio un intersect tra la
                geometria e il suo centroide. Se fossse un'isola
                dovrebbe tornarmi False... o sbaglio?
                
                iface = qgis.utils.iface
                lyr = iface.activeLayer()
                prov = lyr.dataProvider()
                attrlist = prov.attributeIndexes()
                prov.select(attrlist)
                feat = QgsFeature()
                for i in range(lyr.featureCount()):
                    prov.nextFeature(feat)
                    geom = feat.geometry()
                    cent = geom.centroid()
                    dentro = geom.intersects(cent)
                    if (dentro):
                        attributes = feat.attributeMap()
                        print 'La feature %s contiene una geometria
                piena' % attributes[0].toInt()[0]
                    else:
                        print '->>> La feature %s sembra contenere
                un'isola' % attributes[0].toInt()[0]
                
                Giovanni
                        
                Il giorno 15 giugno 2011 11:00, marco zanieri
                <marcozanieri@gmail.com> ha scritto:
                        
                        Salve,
                        ho un problema con alcune geometrie areali,
                        avrei la necessità di verificare il senso di
                        digitalizzazione (orario:interno poligono
                        pieno; antiorario: interno poligono vuoto);
                        esiste quest possibilità in Qgis?
                        
                        Grazie mille,
                        marco
                        
                        --
                                    dott. Marco Zanieri
                           e-mail: marcozanieri@gmail.com
                        
                                   cartografia tematica
                                  banche dati territoriali
                             sistemi informativi geografici
                              applicazioni GIS e webGIS
                        
                        _______________________________________________
                        Iscriviti all'associazione GFOSS.it:
                        http://www.gfoss.it/drupal/iscrizione
                        Gfoss@lists.gfoss.it
                        http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
                        Questa e' una lista di discussione pubblica
                        aperta a tutti.
                        Non inviate messaggi commerciali.
                        I messaggi di questa lista non rispecchiano
                        necessariamente
                        le posizioni dell'Associazione GFOSS.it.
                        518 iscritti al 3.6.2011
                
--
            dott. Marco Zanieri
   e-mail: marcozanieri@gmail.com

           cartografia tematica
          banche dati territoriali
     sistemi informativi geografici
      applicazioni GIS e webGIS

_______________________________________________
Iscriviti all'associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell'Associazione GFOSS.it.
518 iscritti al 3.6.2011

Mah, in PostGIS ho trovato solo:

ST_LineCrossingDirection()

che accetta due linestring come parametro. La seconda va generata dal
congiungendo il punto originale (mypoint) con la proiezione del punto
sulla prima (mypline):

ST_LineCrossingDirection(
  mypline,
  st_makeline(
    mypoint,
    ST_Line_Interpolate_Point(
      mypline,
      ST_Line_Locate_Point(
        mypline,
        mypoint
      )
    )
  )
)

Se restituisce -1, il punto è a sinistra; +1, il punto è a destra.

Mi pare strano non ci sia qualcosa di più alto livello.

Torno al lavoro, la pausa di divertimento è finita :))))

Sig

Il giorno mer, 15/06/2011 alle 14.36 +0200, Luca Sigfrido Percich ha
scritto:

Ciao a tutti,

in teoria si fa così:

- estraggo il linearring del poligono;
- ne genero il centroide
- determino la proiezione e la segmentazione del centroide sulla
polilinea che corrisponde al linearring.

Se il centroide cade in "destra geometrica", allora l'arco è
digitalizzato in senso orario; altrimenti, in senso antiorario.

Per proiezione e segmentazione, intendo ad esempio le funzioni LRS di
PostGIS (ST_Line_Locate_Point() per la segmentazione, ma non trovo
quella che mi dà il lato). Ci sono funzioni analoghe in JTS, quindi in
GEOS e quindi in python. Secoli fa scrissi una classe Java per gestire
la segmentazione usando JTS, se ti serve te la spedisco, ma spero che
qualcuno abbia già implementato una funzione ad alto livello che dato un
punto e una polilinea mi dice a che distanza, da che lato e con che
proiezione questo cade su quella. Le funzioni JTS lavorano a livello di
punto e segmento (2 punti), quindi devi fare un ciclo per lavorare su
una linestring.

Fammi sapere se trovi qualcosa, mi interessa molto (soprattutto mi piace
l'idea di non dover fare il porting di quanto già implementato anni fa e
potermi così dedicare a cose nuove)

Sig

Il giorno mer, 15/06/2011 alle 13.00 +0200, marco zanieri ha scritto:
> ho provato con la tematizzazione delle coordinate e si riesce a
> comprendere l'orientamento...grazie mille Luca; credo comunque che
> possa risultare molto utile la possibilità di individuare la sequenza
> delle coordinate di una geometria e comunque una procedura automatica
> che possa riconoscere la presenza di eventuali isole (pieni e
> vuoti)...
> grazie a tutti
>
> Il giorno 15 giugno 2011 12:49, G. Allegri <giohappy@gmail.com> ha
> scritto:
> Il test del centroide è un test del pirla, scusate, perché il
> centroide può cadere anche fuori dalla geometria...
> Posterò una domanda nell ml di sviluppo ma intanto lo chiedo
> anche qua: come si può sapere l'orientamento di una geometria
> tramite le api di Qgis?
>
> giovanni
>
> Il giorno 15 giugno 2011 12:36, G. Allegri
> <giohappy@gmail.com> ha scritto:
>
>
> Dall'API di QGis non trovo un metodo per determinare
> l'orientamento della sequenza di coordinate... o se si
> tratta di un'isola. Un modo grezzo ma, credo, efficace
> è lanciare questo script nella console di Python
> dentro QGis, nel quale faccio un intersect tra la
> geometria e il suo centroide. Se fossse un'isola
> dovrebbe tornarmi False... o sbaglio?
>
> iface = qgis.utils.iface
> lyr = iface.activeLayer()
> prov = lyr.dataProvider()
> attrlist = prov.attributeIndexes()
> prov.select(attrlist)
> feat = QgsFeature()
> for i in range(lyr.featureCount()):
> prov.nextFeature(feat)
> geom = feat.geometry()
> cent = geom.centroid()
> dentro = geom.intersects(cent)
> if (dentro):
> attributes = feat.attributeMap()
> print 'La feature %s contiene una geometria
> piena' % attributes[0].toInt()[0]
> else:
> print '->>> La feature %s sembra contenere
> un'isola' % attributes[0].toInt()[0]
>
>
> Giovanni
>
>
>
>
>
> Il giorno 15 giugno 2011 11:00, marco zanieri
> <marcozanieri@gmail.com> ha scritto:
>
> Salve,
> ho un problema con alcune geometrie areali,
> avrei la necessità di verificare il senso di
> digitalizzazione (orario:interno poligono
> pieno; antiorario: interno poligono vuoto);
> esiste quest possibilità in Qgis?
>
> Grazie mille,
> marco
>
> --
> dott. Marco Zanieri
> e-mail: marcozanieri@gmail.com
>
> cartografia tematica
> banche dati territoriali
> sistemi informativi geografici
> applicazioni GIS e webGIS
>
>
>
>
>
>
> _______________________________________________
> Iscriviti all'associazione GFOSS.it:
> http://www.gfoss.it/drupal/iscrizione
> Gfoss@lists.gfoss.it
> http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
> Questa e' una lista di discussione pubblica
> aperta a tutti.
> Non inviate messaggi commerciali.
> I messaggi di questa lista non rispecchiano
> necessariamente
> le posizioni dell'Associazione GFOSS.it.
> 518 iscritti al 3.6.2011
>
>
>
>
>
>
>
> --
> dott. Marco Zanieri
> e-mail: marcozanieri@gmail.com
>
> cartografia tematica
> banche dati territoriali
> sistemi informativi geografici
> applicazioni GIS e webGIS
>
>
>
>
> _______________________________________________
> Iscriviti all'associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
> Gfoss@lists.gfoss.it
> http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
> Questa e' una lista di discussione pubblica aperta a tutti.
> Non inviate messaggi commerciali.
> I messaggi di questa lista non rispecchiano necessariamente
> le posizioni dell'Associazione GFOSS.it.
> 518 iscritti al 3.6.2011

_______________________________________________
Iscriviti all'associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell'Associazione GFOSS.it.
518 iscritti al 3.6.2011

Più semplicemte si può usare il metodo dell’area. Se è positiva è antiorario, e viceversa.

Console Python di Qgis:

def senso(p):
area = 0
p0 = p[0]
for p1 in p[1:]:
area += (p0[0] * p1[1]) - (p0[1] * p1[0]);
p0 = p1

if area>0: # antiorario
return 0
else: #orario
return 1

sensi = {0:‘antiorario’,1:‘orario’}
iface = qgis.utils.iface
lyr = iface.activeLayer()
prov = lyr.dataProvider()
attrlist = prov.attributeIndexes()
prov.select(attrlist)
feat = QgsFeature()
for i in range(lyr.featureCount()):
prov.nextFeature(feat)
geom = feat.geometry()
wkbt = geom.wkbType()
attributes = feat.attributeMap()
print ‘->>> Feature GIS: %s’ % attributes[0].toInt()[0]
if wkbt==3:
poly = geom.asPolygon()
#seq_order(poly[0],attributes)
s = senso(poly[0])
print ‘Poligono %s: %s’ % (i,sensi[s])
else:
poly = geom.asMultiPolygon()
for p in poly:
#seq_order(p[0],attributes)
s = senso(p[0])
print ‘Poligono %s: %s’ % (i,sensi[s])

Il giorno 15 giugno 2011 14:58, Luca Sigfrido Percich <sigfrido@tiscali.it> ha scritto:

Mah, in PostGIS ho trovato solo:

ST_LineCrossingDirection()

che accetta due linestring come parametro. La seconda va generata dal
congiungendo il punto originale (mypoint) con la proiezione del punto
sulla prima (mypline):

ST_LineCrossingDirection(
mypline,
st_makeline(
mypoint,
ST_Line_Interpolate_Point(
mypline,
ST_Line_Locate_Point(
mypline,
mypoint
)
)
)
)

Se restituisce -1, il punto è a sinistra; +1, il punto è a destra.

Mi pare strano non ci sia qualcosa di più alto livello.

Torno al lavoro, la pausa di divertimento è finita :))))

Sig

Il giorno mer, 15/06/2011 alle 14.36 +0200, Luca Sigfrido Percich ha
scritto:

Ciao a tutti,

in teoria si fa così:

  • estraggo il linearring del poligono;
  • ne genero il centroide
  • determino la proiezione e la segmentazione del centroide sulla
    polilinea che corrisponde al linearring.

Se il centroide cade in “destra geometrica”, allora l’arco è
digitalizzato in senso orario; altrimenti, in senso antiorario.

Per proiezione e segmentazione, intendo ad esempio le funzioni LRS di
PostGIS (ST_Line_Locate_Point() per la segmentazione, ma non trovo
quella che mi dà il lato). Ci sono funzioni analoghe in JTS, quindi in
GEOS e quindi in python. Secoli fa scrissi una classe Java per gestire
la segmentazione usando JTS, se ti serve te la spedisco, ma spero che
qualcuno abbia già implementato una funzione ad alto livello che dato un
punto e una polilinea mi dice a che distanza, da che lato e con che
proiezione questo cade su quella. Le funzioni JTS lavorano a livello di
punto e segmento (2 punti), quindi devi fare un ciclo per lavorare su
una linestring.

Fammi sapere se trovi qualcosa, mi interessa molto (soprattutto mi piace
l’idea di non dover fare il porting di quanto già implementato anni fa e
potermi così dedicare a cose nuove)

Sig

Il giorno mer, 15/06/2011 alle 13.00 +0200, marco zanieri ha scritto:

ho provato con la tematizzazione delle coordinate e si riesce a
comprendere l’orientamento…grazie mille Luca; credo comunque che
possa risultare molto utile la possibilità di individuare la sequenza
delle coordinate di una geometria e comunque una procedura automatica
che possa riconoscere la presenza di eventuali isole (pieni e
vuoti)…
grazie a tutti

Il giorno 15 giugno 2011 12:49, G. Allegri <giohappy@gmail.com> ha
scritto:
Il test del centroide è un test del pirla, scusate, perché il
centroide può cadere anche fuori dalla geometria…
Posterò una domanda nell ml di sviluppo ma intanto lo chiedo
anche qua: come si può sapere l’orientamento di una geometria
tramite le api di Qgis?

giovanni

Il giorno 15 giugno 2011 12:36, G. Allegri
<giohappy@gmail.com> ha scritto:

Dall’API di QGis non trovo un metodo per determinare
l’orientamento della sequenza di coordinate… o se si
tratta di un’isola. Un modo grezzo ma, credo, efficace
è lanciare questo script nella console di Python
dentro QGis, nel quale faccio un intersect tra la
geometria e il suo centroide. Se fossse un’isola
dovrebbe tornarmi False… o sbaglio?

iface = qgis.utils.iface
lyr = iface.activeLayer()
prov = lyr.dataProvider()
attrlist = prov.attributeIndexes()
prov.select(attrlist)
feat = QgsFeature()
for i in range(lyr.featureCount()):
prov.nextFeature(feat)
geom = feat.geometry()
cent = geom.centroid()
dentro = geom.intersects(cent)
if (dentro):
attributes = feat.attributeMap()
print ‘La feature %s contiene una geometria
piena’ % attributes[0].toInt()[0]
else:
print ‘->>> La feature %s sembra contenere
un’isola’ % attributes[0].toInt()[0]

Giovanni

Il giorno 15 giugno 2011 11:00, marco zanieri
<marcozanieri@gmail.com> ha scritto:

Salve,
ho un problema con alcune geometrie areali,
avrei la necessità di verificare il senso di
digitalizzazione (orario:interno poligono
pieno; antiorario: interno poligono vuoto);
esiste quest possibilità in Qgis?

Grazie mille,
marco


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi geografici
applicazioni GIS e webGIS


Iscriviti all’associazione GFOSS.it:
http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica
aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano
necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi geografici
applicazioni GIS e webGIS


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011

Ciao Giovanni,

grazie, questa dell'area non la sapevo. :wink:

Tra l'altro come notavi giustamente il centroide può cadere fuori se il
poligono è convoluto, a meno che l'API specifica non ti fornisca una
funzione che garantisca la posizione del centroide nel poligono, quindi
il metodo da me proposto è rischioso, si dovrebbe sostituire il
centroide ad un qualunque punto sicuramente interno al poligono.

Sono andato a ripescare una vecchia funzione MapInfo che scrissi per
determinare il verso di rotazione delle coordinate di un linearring:
data la matrice dei punti che lo costituiscono (1..n), e data una
matrice di 4 elementi che contengono gli indici dei punti,
rispettivamente:
1 - y = maxY
2 - x = maxX
3 - y = minY
4 - x = minX

se gli indici contenuti in almeno 3 dei 4 elementi sono crescenti, la
linestring ruota in senso orario.

Forse computazionalmente meno intensivo (nessuna moltiplicazione) ma
decisamente meno elegante :slight_smile:

Proverò a fare il porting in python

Sig

Il giorno mer, 15/06/2011 alle 15.08 +0200, G. Allegri ha scritto:

Più semplicemte si può usare il metodo dell'area. Se è positiva è
antiorario, e viceversa.

Console Python di Qgis:

def senso(p):
    area = 0
    p0 = p[0]
    for p1 in p[1:]:
        area += (p0[0] * p1[1]) - (p0[1] * p1[0]);
        p0 = p1
    
    if area>0: # antiorario
        return 0
    else: #orario
        return 1
    
sensi = {0:'antiorario',1:'orario'}
iface = qgis.utils.iface
lyr = iface.activeLayer()
prov = lyr.dataProvider()
attrlist = prov.attributeIndexes()
prov.select(attrlist)
feat = QgsFeature()
for i in range(lyr.featureCount()):
    prov.nextFeature(feat)
    geom = feat.geometry()
    wkbt = geom.wkbType()
    attributes = feat.attributeMap()
    print '->>> Feature GIS: %s' % attributes[0].toInt()[0]
    if wkbt==3:
        poly = geom.asPolygon()
        #seq_order(poly[0],attributes)
        s = senso(poly[0])
        print 'Poligono %s: %s' % (i,sensi[s])
    else:
        poly = geom.asMultiPolygon()
        for p in poly:
            #seq_order(p[0],attributes)
            s = senso(p[0])
            print 'Poligono %s: %s' % (i,sensi[s])

Il giorno 15 giugno 2011 14:58, Luca Sigfrido Percich
<sigfrido@tiscali.it> ha scritto:
        
        Mah, in PostGIS ho trovato solo:
        
        ST_LineCrossingDirection()
        
        che accetta due linestring come parametro. La seconda va
        generata dal
        congiungendo il punto originale (mypoint) con la proiezione
        del punto
        sulla prima (mypline):
        
        ST_LineCrossingDirection(
               mypline,
               st_makeline(
                       mypoint,
                       ST_Line_Interpolate_Point(
                               mypline,
                               ST_Line_Locate_Point(
                                       mypline,
                                       mypoint
                               )
                       )
               )
        )
        
        Se restituisce -1, il punto è a sinistra; +1, il punto è a
        destra.
        
        Mi pare strano non ci sia qualcosa di più alto livello.
        
        Torno al lavoro, la pausa di divertimento è finita :))))
        
        Sig
        
        Il giorno mer, 15/06/2011 alle 14.36 +0200, Luca Sigfrido
        Percich ha
        scritto:
        
        > Ciao a tutti,
        >
        > in teoria si fa così:
        >
        > - estraggo il linearring del poligono;
        > - ne genero il centroide
        > - determino la proiezione e la segmentazione del centroide
        sulla
        > polilinea che corrisponde al linearring.
        >
        > Se il centroide cade in "destra geometrica", allora l'arco è
        > digitalizzato in senso orario; altrimenti, in senso
        antiorario.
        >
        >
        > Per proiezione e segmentazione, intendo ad esempio le
        funzioni LRS di
        > PostGIS (ST_Line_Locate_Point() per la segmentazione, ma non
        trovo
        > quella che mi dà il lato). Ci sono funzioni analoghe in JTS,
        quindi in
        > GEOS e quindi in python. Secoli fa scrissi una classe Java
        per gestire
        > la segmentazione usando JTS, se ti serve te la spedisco, ma
        spero che
        > qualcuno abbia già implementato una funzione ad alto livello
        che dato un
        > punto e una polilinea mi dice a che distanza, da che lato e
        con che
        > proiezione questo cade su quella. Le funzioni JTS lavorano a
        livello di
        > punto e segmento (2 punti), quindi devi fare un ciclo per
        lavorare su
        > una linestring.
        >
        > Fammi sapere se trovi qualcosa, mi interessa molto
        (soprattutto mi piace
        > l'idea di non dover fare il porting di quanto già
        implementato anni fa e
        > potermi così dedicare a cose nuove)
        >
        > Sig
        >
        >
        > Il giorno mer, 15/06/2011 alle 13.00 +0200, marco zanieri ha
        scritto:
        > > ho provato con la tematizzazione delle coordinate e si
        riesce a
        > > comprendere l'orientamento...grazie mille Luca; credo
        comunque che
        > > possa risultare molto utile la possibilità di individuare
        la sequenza
        > > delle coordinate di una geometria e comunque una procedura
        automatica
        > > che possa riconoscere la presenza di eventuali isole
        (pieni e
        > > vuoti)...
        > > grazie a tutti
        > >
        > > Il giorno 15 giugno 2011 12:49, G. Allegri
        <giohappy@gmail.com> ha
        > > scritto:
        > > Il test del centroide è un test del pirla,
        scusate, perché il
        > > centroide può cadere anche fuori dalla
        geometria...
        > > Posterò una domanda nell ml di sviluppo ma intanto
        lo chiedo
        > > anche qua: come si può sapere l'orientamento di
        una geometria
        > > tramite le api di Qgis?
        > >
        > > giovanni
        > >
        > > Il giorno 15 giugno 2011 12:36, G. Allegri
        > > <giohappy@gmail.com> ha scritto:
        > >
        > >
        > > Dall'API di QGis non trovo un metodo per
        determinare
        > > l'orientamento della sequenza di
        coordinate... o se si
        > > tratta di un'isola. Un modo grezzo ma,
        credo, efficace
        > > è lanciare questo script nella console di
        Python
        > > dentro QGis, nel quale faccio un intersect
        tra la
        > > geometria e il suo centroide. Se fossse
        un'isola
        > > dovrebbe tornarmi False... o sbaglio?
        > >
        > > iface = qgis.utils.iface
        > > lyr = iface.activeLayer()
        > > prov = lyr.dataProvider()
        > > attrlist = prov.attributeIndexes()
        > > prov.select(attrlist)
        > > feat = QgsFeature()
        > > for i in range(lyr.featureCount()):
        > > prov.nextFeature(feat)
        > > geom = feat.geometry()
        > > cent = geom.centroid()
        > > dentro = geom.intersects(cent)
        > > if (dentro):
        > > attributes = feat.attributeMap()
        > > print 'La feature %s contiene una
        geometria
        > > piena' % attributes[0].toInt()[0]
        > > else:
        > > print '->>> La feature %s sembra
        contenere
        > > un'isola' % attributes[0].toInt()[0]
        > >
        > >
        > > Giovanni
        > >
        > >
        > >
        > >
        > >
        > > Il giorno 15 giugno 2011 11:00, marco
        zanieri
        > > <marcozanieri@gmail.com> ha scritto:
        > >
        > > Salve,
        > > ho un problema con alcune
        geometrie areali,
        > > avrei la necessità di verificare
        il senso di
        > > digitalizzazione (orario:interno
        poligono
        > > pieno; antiorario: interno
        poligono vuoto);
        > > esiste quest possibilità in Qgis?
        > >
        > > Grazie mille,
        > > marco
        > >
        > > --
        > > dott. Marco Zanieri
        > > e-mail: marcozanieri@gmail.com
        > >
        > > cartografia tematica
        > > banche dati territoriali
        > > sistemi informativi
        geografici
        > > applicazioni GIS e webGIS
        > >
        > >
        > >
        > >
        > >
        > >
        > >
        _______________________________________________
        > > Iscriviti all'associazione
        GFOSS.it:
        > >
        http://www.gfoss.it/drupal/iscrizione
        > > Gfoss@lists.gfoss.it
        > >
        http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
        > > Questa e' una lista di discussione
        pubblica
        > > aperta a tutti.
        > > Non inviate messaggi commerciali.
        > > I messaggi di questa lista non
        rispecchiano
        > > necessariamente
        > > le posizioni dell'Associazione
        GFOSS.it.
        > > 518 iscritti al 3.6.2011
        > >
        > >
        > >
        > >
        > >
        > >
        > >
        > > --
        > > dott. Marco Zanieri
        > > e-mail: marcozanieri@gmail.com
        > >
        > > cartografia tematica
        > > banche dati territoriali
        > > sistemi informativi geografici
        > > applicazioni GIS e webGIS
        > >
        > >
        > >
        > >
        > > _______________________________________________
        > > Iscriviti all'associazione GFOSS.it:
        http://www.gfoss.it/drupal/iscrizione
        > > Gfoss@lists.gfoss.it
        > > http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
        > > Questa e' una lista di discussione pubblica aperta a
        tutti.
        > > Non inviate messaggi commerciali.
        > > I messaggi di questa lista non rispecchiano
        necessariamente
        > > le posizioni dell'Associazione GFOSS.it.
        > > 518 iscritti al 3.6.2011
        >
        > _______________________________________________
        > Iscriviti all'associazione GFOSS.it:
        http://www.gfoss.it/drupal/iscrizione
        > Gfoss@lists.gfoss.it
        > http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
        > Questa e' una lista di discussione pubblica aperta a tutti.
        > Non inviate messaggi commerciali.
        > I messaggi di questa lista non rispecchiano necessariamente
        > le posizioni dell'Associazione GFOSS.it.
        > 518 iscritti al 3.6.2011
        
        _______________________________________________
        Iscriviti all'associazione GFOSS.it:
        http://www.gfoss.it/drupal/iscrizione
        Gfoss@lists.gfoss.it
        http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
        Questa e' una lista di discussione pubblica aperta a tutti.
        Non inviate messaggi commerciali.
        I messaggi di questa lista non rispecchiano necessariamente
        le posizioni dell'Associazione GFOSS.it.
        518 iscritti al 3.6.2011

Il giorno 15 giugno 2011 15:21, Luca Sigfrido Percich <sigfrido@tiscali.it> ha scritto:

Ciao Giovanni,

grazie, questa dell’area non la sapevo. :wink:

Tra l’altro come notavi giustamente il centroide può cadere fuori se il
poligono è convoluto, a meno che l’API specifica non ti fornisca una
funzione che garantisca la posizione del centroide nel poligono, quindi
il metodo da me proposto è rischioso, si dovrebbe sostituire il
centroide ad un qualunque punto sicuramente interno al poligono.

Sono andato a ripescare una vecchia funzione MapInfo che scrissi per
determinare il verso di rotazione delle coordinate di un linearring:
data la matrice dei punti che lo costituiscono (1…n), e data una
matrice di 4 elementi che contengono gli indici dei punti,
rispettivamente:
1 - y = maxY
2 - x = maxX
3 - y = minY
4 - x = minX

non ho capito cosa contiene la matrice dei 4 indici…

se gli indici contenuti in almeno 3 dei 4 elementi sono crescenti, la
linestring ruota in senso orario.

Forse computazionalmente meno intensivo (nessuna moltiplicazione) ma
decisamente meno elegante :slight_smile:

Proverò a fare il porting in python

Sig

Il giorno mer, 15/06/2011 alle 15.08 +0200, G. Allegri ha scritto:

Più semplicemte si può usare il metodo dell’area. Se è positiva è
antiorario, e viceversa.

Console Python di Qgis:

def senso(p):
area = 0
p0 = p[0]
for p1 in p[1:]:
area += (p0[0] * p1[1]) - (p0[1] * p1[0]);
p0 = p1

if area>0: # antiorario
return 0
else: #orario
return 1

sensi = {0:‘antiorario’,1:‘orario’}
iface = qgis.utils.iface
lyr = iface.activeLayer()
prov = lyr.dataProvider()
attrlist = prov.attributeIndexes()
prov.select(attrlist)
feat = QgsFeature()
for i in range(lyr.featureCount()):
prov.nextFeature(feat)
geom = feat.geometry()
wkbt = geom.wkbType()
attributes = feat.attributeMap()
print ‘->>> Feature GIS: %s’ % attributes[0].toInt()[0]
if wkbt==3:
poly = geom.asPolygon()
#seq_order(poly[0],attributes)
s = senso(poly[0])
print ‘Poligono %s: %s’ % (i,sensi[s])
else:
poly = geom.asMultiPolygon()
for p in poly:
#seq_order(p[0],attributes)
s = senso(p[0])
print ‘Poligono %s: %s’ % (i,sensi[s])

Il giorno 15 giugno 2011 14:58, Luca Sigfrido Percich
<sigfrido@tiscali.it> ha scritto:

Mah, in PostGIS ho trovato solo:

ST_LineCrossingDirection()

che accetta due linestring come parametro. La seconda va
generata dal
congiungendo il punto originale (mypoint) con la proiezione
del punto
sulla prima (mypline):

ST_LineCrossingDirection(
mypline,
st_makeline(
mypoint,
ST_Line_Interpolate_Point(
mypline,
ST_Line_Locate_Point(
mypline,
mypoint
)
)
)
)

Se restituisce -1, il punto è a sinistra; +1, il punto è a
destra.

Mi pare strano non ci sia qualcosa di più alto livello.

Torno al lavoro, la pausa di divertimento è finita :))))

Sig

Il giorno mer, 15/06/2011 alle 14.36 +0200, Luca Sigfrido
Percich ha
scritto:

Ciao a tutti,

in teoria si fa così:

  • estraggo il linearring del poligono;
  • ne genero il centroide
  • determino la proiezione e la segmentazione del centroide
    sulla
    polilinea che corrisponde al linearring.

Se il centroide cade in “destra geometrica”, allora l’arco è
digitalizzato in senso orario; altrimenti, in senso
antiorario.

Per proiezione e segmentazione, intendo ad esempio le
funzioni LRS di
PostGIS (ST_Line_Locate_Point() per la segmentazione, ma non
trovo
quella che mi dà il lato). Ci sono funzioni analoghe in JTS,
quindi in
GEOS e quindi in python. Secoli fa scrissi una classe Java
per gestire
la segmentazione usando JTS, se ti serve te la spedisco, ma
spero che
qualcuno abbia già implementato una funzione ad alto livello
che dato un
punto e una polilinea mi dice a che distanza, da che lato e
con che
proiezione questo cade su quella. Le funzioni JTS lavorano a
livello di
punto e segmento (2 punti), quindi devi fare un ciclo per
lavorare su
una linestring.

Fammi sapere se trovi qualcosa, mi interessa molto
(soprattutto mi piace
l’idea di non dover fare il porting di quanto già
implementato anni fa e
potermi così dedicare a cose nuove)

Sig

Il giorno mer, 15/06/2011 alle 13.00 +0200, marco zanieri ha
scritto:

ho provato con la tematizzazione delle coordinate e si
riesce a
comprendere l’orientamento…grazie mille Luca; credo
comunque che
possa risultare molto utile la possibilità di individuare
la sequenza
delle coordinate di una geometria e comunque una procedura
automatica
che possa riconoscere la presenza di eventuali isole
(pieni e
vuoti)…
grazie a tutti

Il giorno 15 giugno 2011 12:49, G. Allegri
<giohappy@gmail.com> ha
scritto:
Il test del centroide è un test del pirla,
scusate, perché il
centroide può cadere anche fuori dalla
geometria…
Posterò una domanda nell ml di sviluppo ma intanto
lo chiedo
anche qua: come si può sapere l’orientamento di
una geometria
tramite le api di Qgis?

giovanni

Il giorno 15 giugno 2011 12:36, G. Allegri
<giohappy@gmail.com> ha scritto:

Dall’API di QGis non trovo un metodo per
determinare
l’orientamento della sequenza di
coordinate… o se si
tratta di un’isola. Un modo grezzo ma,
credo, efficace
è lanciare questo script nella console di
Python
dentro QGis, nel quale faccio un intersect
tra la
geometria e il suo centroide. Se fossse
un’isola
dovrebbe tornarmi False… o sbaglio?

iface = qgis.utils.iface
lyr = iface.activeLayer()
prov = lyr.dataProvider()
attrlist = prov.attributeIndexes()
prov.select(attrlist)
feat = QgsFeature()
for i in range(lyr.featureCount()):
prov.nextFeature(feat)
geom = feat.geometry()
cent = geom.centroid()
dentro = geom.intersects(cent)
if (dentro):
attributes = feat.attributeMap()
print ‘La feature %s contiene una
geometria
piena’ % attributes[0].toInt()[0]
else:
print ‘->>> La feature %s sembra
contenere
un’isola’ % attributes[0].toInt()[0]

Giovanni

Il giorno 15 giugno 2011 11:00, marco
zanieri
<marcozanieri@gmail.com> ha scritto:

Salve,
ho un problema con alcune
geometrie areali,
avrei la necessità di verificare
il senso di
digitalizzazione (orario:interno
poligono
pieno; antiorario: interno
poligono vuoto);
esiste quest possibilità in Qgis?

Grazie mille,
marco


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi
geografici
applicazioni GIS e webGIS


Iscriviti all’associazione
GFOSS.it:

http://www.gfoss.it/drupal/iscrizione

Gfoss@lists.gfoss.it

http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss

Questa e’ una lista di discussione
pubblica
aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non
rispecchiano
necessariamente
le posizioni dell’Associazione
GFOSS.it.
518 iscritti al 3.6.2011


dott. Marco Zanieri
e-mail: marcozanieri@gmail.com

cartografia tematica
banche dati territoriali
sistemi informativi geografici
applicazioni GIS e webGIS


Iscriviti all’associazione GFOSS.it:
http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a
tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano
necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011


Iscriviti all’associazione GFOSS.it:
http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011


Iscriviti all’associazione GFOSS.it:
http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011

Scusa, lo so, mi sono espresso frettolosamente.

Se p[0]..p[n] è la matrice dei punti che compongono il poligono ordinata
secondo il verso geometrico,

bound[0]..bound[3] è una matrice che contiene 4 indici numerici tali
che:

p[bound[0]] è il punto a ordinata massima (p[bound[0]].y = y_min)
p[bound[1]] è il punto ad ascissa massima
p[bound[2]] è il punto as ordinata minima
p[bound[3]] è il punto ad ascissa minima

bound va determinata con una scansione di p da 0 a n (n - 1 se il
primo e l'ultimo punto del linearring coincidono), confrontando ascisse
e ordinate dei singoli punti con quelle di due punti max.x, max.y,
min.x, min.y.

I p[bound] sono i punti in cui il poligono tocca il bounding box.
L'algoritmo funziona anche se ho più punti che intersecano lo steso lato
del bbox, in questo caso viene sempre preso l'ultimo trovato.

Dati dati tre indici a, b, c con
  a >= 0
  c > b > a
  c < 4

se bound[a] < bound[b] < bound[c]

allora il ring è digitalizzato in senso orario.

Scusami ma sono appena agli inizi con python, non so fornirti un esempio
sintatticamente corretto.

Sig

Il giorno mer, 15/06/2011 alle 15.24 +0200, G. Allegri ha scritto:

non ho capito cosa contiene la matrice dei 4 indici....

Interessante, questo metodo non lo conoscevo, grazie :wink:

Il giorno 15 giugno 2011 15:46, Luca Sigfrido Percich <sigfrido@tiscali.it> ha scritto:

Scusa, lo so, mi sono espresso frettolosamente.

Se p[0]…p[n] è la matrice dei punti che compongono il poligono ordinata
secondo il verso geometrico,

bound[0]…bound[3] è una matrice che contiene 4 indici numerici tali
che:

p[bound[0]] è il punto a ordinata massima (p[bound[0]].y = y_min)
p[bound[1]] è il punto ad ascissa massima
p[bound[2]] è il punto as ordinata minima
p[bound[3]] è il punto ad ascissa minima

bound va determinata con una scansione di p da 0 a n (n - 1 se il
primo e l’ultimo punto del linearring coincidono), confrontando ascisse
e ordinate dei singoli punti con quelle di due punti max.x, max.y,
min.x, min.y.

I p[bound] sono i punti in cui il poligono tocca il bounding box.
L’algoritmo funziona anche se ho più punti che intersecano lo steso lato
del bbox, in questo caso viene sempre preso l’ultimo trovato.

Dati dati tre indici a, b, c con
a >= 0
c > b > a
c < 4

se bound[a] < bound[b] < bound[c]

allora il ring è digitalizzato in senso orario.

Scusami ma sono appena agli inizi con python, non so fornirti un esempio
sintatticamente corretto.

Sig

Il giorno mer, 15/06/2011 alle 15.24 +0200, G. Allegri ha scritto:

non ho capito cosa contiene la matrice dei 4 indici…

...che dici, apriamo un giro di scommesse sui risultati dei benchmark?

Gfoss: Gambling on free open source software

:)))))))

Sig

PS: prima devo imparare python, però. :o|

Il giorno mer, 15/06/2011 alle 15.54 +0200, G. Allegri ha scritto:

Interessante, questo metodo non lo conoscevo, grazie :wink:

On Wed, 15 Jun 2011 15:46:29 +0200
Luca Sigfrido Percich <sigfrido@tiscali.it> wrote:

Scusa, lo so, mi sono espresso frettolosamente.

Se p[0]..p[n] è la matrice dei punti che compongono il poligono
ordinata secondo il verso geometrico,

bound[0]..bound[3] è una matrice che contiene 4 indici numerici tali
che:

p[bound[0]] è il punto a ordinata massima (p[bound[0]].y = y_min)
p[bound[1]] è il punto ad ascissa massima
p[bound[2]] è il punto as ordinata minima
p[bound[3]] è il punto ad ascissa minima
........

ciao Luca,

0) mi sembra che il tuo algoritmo si basi sul fatto che per passare con
un orientamento ORARIO dall'ordinata max a quella min occorre passare
prima per l'ascissa max (o per passare dall'ascissa max a quella min
occorre prima passare dall'ordinata min): questo e' il senso
dell'ordine crescente richiesto per i dati secondo i tre indici 0,1,2 o
1,2,3;

1) mi sembra che l'algoritmo quindi funzioni solo per poligoni non
intrecciati (non ho un background geografico - cartografico, per cui ne
parlo in termini puramente geometrici);

2) non e' robusto: mi sembra che nel caso di un poligono antiorario
(1,1), (3,2), (4,4), (2,3), (1,1) o orario (1,1), (2,3), (4,4), (3,2),
(1,1) non sia in grado di determinarne il senso;

3) curiosita': quale la fonte di questo algoritmo?

4) il metodo dell'area proposto da Giovanni e' basato
sulla somma algebrica delle aree sottese ad ogni lato (cosi' mi sembra
di ricordare nelle versioni che ho visto pubblicate: Newmann/Sproull?
Foley/VanDam? Knuth? se e' importante cerco di trovarlo) quindi dovrebbe
dare comunque il valore dell'area; il segno sara' positivo o negativo a
seconda del senso di percorrenza (dovrebbe funzionare anche
per poligoni intrecciati);

5) un concetto analogo e' quello del momento statico di una forza
(un lato) rispetto ad un punto che puo' destrogiro o sinistrogiro, ma
occorre verificare e inoltre il calcolo probabilmente e' molto vicino a
quello fatto con il metodo dell'area;

ciao,
giuliano

Ciao a tutti,

ciao Giuliano e grazie per esserti preso la briga di verificare gli
algoritmi!

Il giorno gio, 16/06/2011 alle 09.48 +0200, giuliano ha scritto:

On Wed, 15 Jun 2011 15:46:29 +0200

0) mi sembra che il tuo algoritmo si basi sul fatto che per passare con
un orientamento ORARIO dall'ordinata max a quella min occorre passare
prima per l'ascissa max (o per passare dall'ascissa max a quella min
occorre prima passare dall'ordinata min): questo e' il senso
dell'ordine crescente richiesto per i dati secondo i tre indici 0,1,2 o
1,2,3;

Esattamente

1) mi sembra che l'algoritmo quindi funzioni solo per poligoni non
intrecciati (non ho un background geografico - cartografico, per cui ne
parlo in termini puramente geometrici);

Esatto; un poligono come quello da te descritto è topologicamente
sbagliato (presenta una autointersezione). Quindi prima di applicare
l'algoritmo, la geometria va controllata. Del resto, un poligono a forma
di 8 in che senso gira? Un anello in orario, e l'altro in antiorario, ma
per il poligono in se il senso di rotazione non ha valore.

A meno che non si desideri conoscere il senso di rotazione in
corrispondenza di un punto esterno al perimetro, ossia il verso del
vettore tangente al poligono nella proiezione del punto sul perimetro.
Questo potrebbe avere senso, ed è anche semplice da implementare (a
occhio).

2) non e' robusto: mi sembra che nel caso di un poligono antiorario
(1,1), (3,2), (4,4), (2,3), (1,1) o orario (1,1), (2,3), (4,4), (3,2),
(1,1) non sia in grado di determinarne il senso;

Hai ragione, hai trovato il caso limite, ovvero quello in due soli punti
soddisfano le 4 condizioni, ovvero quando due punti coincidono con 2
vertici opposti dell'MBR. Ho verificato il mio codice originale e manca
il controllo, che potrebbe essere implementato così (ora non ho molto
tempo per ragionarci bene):

Se al termine del ciclo sui vertici, la matrice bound contiene solo 2
indici che si ripetono, quindi nel tuo esempio

p[bound[0]] = 3 (4,4) ordinata max
p[bound[1]] = 3 (4, 4) ascissa max
p[bound[2]] = 1 (1, 1) ordinata minima
p[bound[3]] = 1 (1, 1) ascissa minima

dovrebbe essere sufficiente cercare un altro vertice per una qualunque
delle condizioni, escludendo dal suo proprio criterio i punti i cui
indici siano già contenuti nella matrice, ovvero:

p[bound[0]] = punto a ordinata massima che non sia 3 o 1 = 4 (2, 3)

Il che risulta in un ordine antiorario (4 3 3 1)

Se scelgo una delle altre condizioni:

p[bound[1]] = 2 (3 2 1 1) ascissa max <> 3 OK
p[bound[2]] = 2 (3 3 2 1) ordinata minima <> 1 OK
p[bound[3]] = 4 (3 3 1 4) ascissa minima <> 1 ?!?!?!

La quarta condizione è problematica: non so se basti "ruotare la
sequenza", ovvero se 3314 = 4331, oppure se la scelta di quale punto
cercare non possa essere casuale ma dipenda da determinate condizioni.
Ci penserò! Dammi una mano se vuoi/puoi.

3) curiosita': quale la fonte di questo algoritmo?

L'ho scritto io, in MapBasic, nel 2006.

4) il metodo dell'area proposto da Giovanni e' basato
sulla somma algebrica delle aree sottese ad ogni lato (cosi' mi sembra
di ricordare nelle versioni che ho visto pubblicate: Newmann/Sproull?
Foley/VanDam? Knuth? se e' importante cerco di trovarlo) quindi dovrebbe
dare comunque il valore dell'area; il segno sara' positivo o negativo a
seconda del senso di percorrenza (dovrebbe funzionare anche
per poligoni intrecciati);

5) un concetto analogo e' quello del momento statico di una forza
(un lato) rispetto ad un punto che puo' destrogiro o sinistrogiro, ma
occorre verificare e inoltre il calcolo probabilmente e' molto vicino a
quello fatto con il metodo dell'area;

Grazie per la spiegazione! Qui fatico un po', purtroppo il mio
background accademico da agronomo non mi aiuta molto.

Vorrei implementarlo in python sulla falsariga dell'implementazione
proposta da Giovanni, ovviamente appena riesco a prendermi il tempo per
imparare python. Nel frattempo, se a qualcuno interessa posso postare il
codice MapBasic.

Ciao e grazie ancora

Sig

ciao,
giuliano

_______________________________________________
Iscriviti all'associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell'Associazione GFOSS.it.
518 iscritti al 3.6.2011

Il risultato interno del metodo dell’area che ho postato ovviamente dovrebbe essere diviso per 2 per avere l’area, ma a noi interessa il segno, quindi evitiamo una divisione inutile.
Qui c’è una delle tante descrizioni [1]

Il discorso sul versore del braccio potrebbe essere implementato calcolando il prodotto vettoriale tra due edge consecutivi, e valutarne poi il segno.

giovanni

[1] http://paulbourke.net/geometry/polyarea/

Il giorno 16 giugno 2011 12:54, Luca Sigfrido Percich <sigfrido@tiscali.it> ha scritto:

Ciao a tutti,

ciao Giuliano e grazie per esserti preso la briga di verificare gli
algoritmi!

Il giorno gio, 16/06/2011 alle 09.48 +0200, giuliano ha scritto:

On Wed, 15 Jun 2011 15:46:29 +0200

  1. mi sembra che il tuo algoritmo si basi sul fatto che per passare con
    un orientamento ORARIO dall’ordinata max a quella min occorre passare
    prima per l’ascissa max (o per passare dall’ascissa max a quella min
    occorre prima passare dall’ordinata min): questo e’ il senso
    dell’ordine crescente richiesto per i dati secondo i tre indici 0,1,2 o
    1,2,3;

Esattamente

  1. mi sembra che l’algoritmo quindi funzioni solo per poligoni non
    intrecciati (non ho un background geografico - cartografico, per cui ne
    parlo in termini puramente geometrici);

Esatto; un poligono come quello da te descritto è topologicamente
sbagliato (presenta una autointersezione). Quindi prima di applicare
l’algoritmo, la geometria va controllata. Del resto, un poligono a forma
di 8 in che senso gira? Un anello in orario, e l’altro in antiorario, ma
per il poligono in se il senso di rotazione non ha valore.

A meno che non si desideri conoscere il senso di rotazione in
corrispondenza di un punto esterno al perimetro, ossia il verso del
vettore tangente al poligono nella proiezione del punto sul perimetro.
Questo potrebbe avere senso, ed è anche semplice da implementare (a
occhio).

  1. non e’ robusto: mi sembra che nel caso di un poligono antiorario
    (1,1), (3,2), (4,4), (2,3), (1,1) o orario (1,1), (2,3), (4,4), (3,2),
    (1,1) non sia in grado di determinarne il senso;

Hai ragione, hai trovato il caso limite, ovvero quello in due soli punti
soddisfano le 4 condizioni, ovvero quando due punti coincidono con 2
vertici opposti dell’MBR. Ho verificato il mio codice originale e manca
il controllo, che potrebbe essere implementato così (ora non ho molto
tempo per ragionarci bene):

Se al termine del ciclo sui vertici, la matrice bound contiene solo 2
indici che si ripetono, quindi nel tuo esempio

p[bound[0]] = 3 (4,4) ordinata max
p[bound[1]] = 3 (4, 4) ascissa max
p[bound[2]] = 1 (1, 1) ordinata minima
p[bound[3]] = 1 (1, 1) ascissa minima

dovrebbe essere sufficiente cercare un altro vertice per una qualunque
delle condizioni, escludendo dal suo proprio criterio i punti i cui
indici siano già contenuti nella matrice, ovvero:

p[bound[0]] = punto a ordinata massima che non sia 3 o 1 = 4 (2, 3)

Il che risulta in un ordine antiorario (4 3 3 1)

Se scelgo una delle altre condizioni:

p[bound[1]] = 2 (3 2 1 1) ascissa max <> 3 OK
p[bound[2]] = 2 (3 3 2 1) ordinata minima <> 1 OK
p[bound[3]] = 4 (3 3 1 4) ascissa minima <> 1 ?!?!?!

La quarta condizione è problematica: non so se basti “ruotare la
sequenza”, ovvero se 3314 = 4331, oppure se la scelta di quale punto
cercare non possa essere casuale ma dipenda da determinate condizioni.
Ci penserò! Dammi una mano se vuoi/puoi.

  1. curiosita’: quale la fonte di questo algoritmo?

L’ho scritto io, in MapBasic, nel 2006.

  1. il metodo dell’area proposto da Giovanni e’ basato
    sulla somma algebrica delle aree sottese ad ogni lato (cosi’ mi sembra
    di ricordare nelle versioni che ho visto pubblicate: Newmann/Sproull?
    Foley/VanDam? Knuth? se e’ importante cerco di trovarlo) quindi dovrebbe
    dare comunque il valore dell’area; il segno sara’ positivo o negativo a
    seconda del senso di percorrenza (dovrebbe funzionare anche
    per poligoni intrecciati);

  2. un concetto analogo e’ quello del momento statico di una forza
    (un lato) rispetto ad un punto che puo’ destrogiro o sinistrogiro, ma
    occorre verificare e inoltre il calcolo probabilmente e’ molto vicino a
    quello fatto con il metodo dell’area;

Grazie per la spiegazione! Qui fatico un po’, purtroppo il mio
background accademico da agronomo non mi aiuta molto.

Vorrei implementarlo in python sulla falsariga dell’implementazione
proposta da Giovanni, ovviamente appena riesco a prendermi il tempo per
imparare python. Nel frattempo, se a qualcuno interessa posso postare il
codice MapBasic.

Ciao e grazie ancora

Sig

ciao,
giuliano


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011

On Thu, 16 Jun 2011 12:54:55 +0200
Luca Sigfrido Percich <sigfrido@tiscali.it> wrote:

Ciao a tutti,

ciao Giuliano e grazie per esserti preso la briga di verificare gli
algoritmi!

mi piacerebbe contribuire :-)))

......

> 2) non e' robusto: mi sembra che nel caso di un poligono antiorario
> (1,1), (3,2), (4,4), (2,3), (1,1) o orario (1,1), (2,3), (4,4),
> (3,2), (1,1) non sia in grado di determinarne il senso;

Hai ragione, hai trovato il caso limite, ovvero quello in due soli
punti soddisfano le 4 condizioni, ovvero quando due punti coincidono
con 2 vertici opposti dell'MBR. Ho verificato il mio codice originale
e manca il controllo, che potrebbe essere implementato così (ora non
ho molto tempo per ragionarci bene):

......

dovrebbe essere sufficiente cercare un altro vertice per una qualunque
delle condizioni, escludendo dal suo proprio criterio i punti i cui
indici siano già contenuti nella matrice, ovvero:

p[bound[0]] = punto a ordinata massima che non sia 3 o 1 = 4 (2, 3)

Il che risulta in un ordine antiorario (4 3 3 1)

.....

temo non basti: se avessi (1,1), (3,3.1), (4,4), (2,3), (1,1)
troveresti il nodo 2 al posto del 4 e riceveresti una informazione
errata:-(( (scusa, da verificare bene, l'ho buttata li' :-))

Ci penserò! Dammi una mano se vuoi/puoi.

con infinito piacere se posso;

> ......

Vorrei implementarlo in python sulla falsariga dell'implementazione
proposta da Giovanni, ovviamente appena riesco a prendermi il tempo
per imparare python. Nel frattempo, se a qualcuno interessa posso
postare il codice MapBasic.

anch'io conosco pochissimo python e nulla della struttura dati di qGIS;
se qualcuno ritenesse non buttato il tempo per qualche dritta diretta
o qualche rinvio a documentazione (nella speranza che questa non
richieda di essere gia' dei guru per capirla come ahime' alcune volte
succede....:slight_smile: si puo' provare;

On Thu, 16 Jun 2011 13:03:16 +0200
"G. Allegri" <giohappy@gmail.com> wrote:

.......

[1] http://paulbourke.net/geometry/polyarea/

visto, grazie;

ciao,
giuliano

Io ho esperienza di sviluppo su Qgis, ma ora non ho proprio tempo libero…
Questo comunque potrebbe esservi utile: http://www.qgis.org/pyqgis-cookbook/

in ogni caso sono a disposizione per darvi una mano, tra un caffè e l’altro,
giovanni

Il giorno 16 giugno 2011 18:42, giuliano <giulianc@tiscali.it> ha scritto:

On Thu, 16 Jun 2011 12:54:55 +0200

Luca Sigfrido Percich <sigfrido@tiscali.it> wrote:

Ciao a tutti,

ciao Giuliano e grazie per esserti preso la briga di verificare gli
algoritmi!

mi piacerebbe contribuire :-)))

  1. non e’ robusto: mi sembra che nel caso di un poligono antiorario
    (1,1), (3,2), (4,4), (2,3), (1,1) o orario (1,1), (2,3), (4,4),
    (3,2), (1,1) non sia in grado di determinarne il senso;

Hai ragione, hai trovato il caso limite, ovvero quello in due soli
punti soddisfano le 4 condizioni, ovvero quando due punti coincidono
con 2 vertici opposti dell’MBR. Ho verificato il mio codice originale
e manca il controllo, che potrebbe essere implementato così (ora non
ho molto tempo per ragionarci bene):

dovrebbe essere sufficiente cercare un altro vertice per una qualunque
delle condizioni, escludendo dal suo proprio criterio i punti i cui
indici siano già contenuti nella matrice, ovvero:

p[bound[0]] = punto a ordinata massima che non sia 3 o 1 = 4 (2, 3)

Il che risulta in un ordine antiorario (4 3 3 1)

temo non basti: se avessi (1,1), (3,3.1), (4,4), (2,3), (1,1)
troveresti il nodo 2 al posto del 4 e riceveresti una informazione
errata:-(( (scusa, da verificare bene, l’ho buttata li’ :-))

Ci penserò! Dammi una mano se vuoi/puoi.

con infinito piacere se posso;

Vorrei implementarlo in python sulla falsariga dell’implementazione
proposta da Giovanni, ovviamente appena riesco a prendermi il tempo
per imparare python. Nel frattempo, se a qualcuno interessa posso
postare il codice MapBasic.

anch’io conosco pochissimo python e nulla della struttura dati di qGIS;
se qualcuno ritenesse non buttato il tempo per qualche dritta diretta
o qualche rinvio a documentazione (nella speranza che questa non
richieda di essere gia’ dei guru per capirla come ahime’ alcune volte
succede…:slight_smile: si puo’ provare;

On Thu, 16 Jun 2011 13:03:16 +0200
“G. Allegri” <giohappy@gmail.com> wrote:

[1] http://paulbourke.net/geometry/polyarea/

visto, grazie;

ciao,
giuliano


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
518 iscritti al 3.6.2011

Il giorno gio, 16/06/2011 alle 18.42 +0200, giuliano ha scritto:

> ......
>
> dovrebbe essere sufficiente cercare un altro vertice per una qualunque
> delle condizioni, escludendo dal suo proprio criterio i punti i cui
> indici siano già contenuti nella matrice, ovvero:
>
> p[bound[0]] = punto a ordinata massima che non sia 3 o 1 = 4 (2, 3)
>
> Il che risulta in un ordine antiorario (4 3 3 1)
>
> .....

temo non basti: se avessi (1,1), (3,3.1), (4,4), (2,3), (1,1)
troveresti il nodo 2 al posto del 4 e riceveresti una informazione
errata:-(( (scusa, da verificare bene, l'ho buttata li' :-))

Hai ragione, la soluzione dev'essere nella scelta di quali elementi di
bounds sostituire.

Oppure vietiamo l'uso dei poligoni con vertici coincidenti con 2 angoli
opposti del proprio MBR. :slight_smile:

anch'io conosco pochissimo python e nulla della struttura dati di qGIS;
se qualcuno ritenesse non buttato il tempo per qualche dritta diretta
o qualche rinvio a documentazione (nella speranza che questa non
richieda di essere gia' dei guru per capirla come ahime' alcune volte
succede....:slight_smile: si puo' provare;

Mi sembra che ci sia molta documentazione, poi c'è PluginBuilder (un
plugin python che può generare lo scheletro di un plugin python), e ci
sono tantissimi plugins tra cui curiosare volendo imparare da esempi
concreti.

Ciao

Sig