[Gfoss] verifica senso digitalizzazione geometria
Luca Sigfrido Percich
sigfrido a tiscali.it
Mer 15 Giu 2011 15:21:35 CEST
Ciao Giovanni,
grazie, questa dell'area non la sapevo. ;)
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 :)
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 a 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 a 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 a 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 a 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 a 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 a 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 a 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 a 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 a 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 a 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
>
>
Maggiori informazioni sulla lista
Gfoss