7 votos

Punto-en-polígono usando Python y OGR

Quiero usar OGR en Python para escribir una simple prueba de punto en polígono:

def point_in_polygon(point, polygon):
    """
    point : [longitude, latitude]

    polygon : [(lon1, lat1), (lon2, lat2), ..., (lonn, latn)]

    """
    # Create spatialReference
    spatialReference = osgeo.osr.SpatialReference()
    spatialReference.SetWellKnownGeogCS("WGS84")
    # Create ring
    ring = osgeo.ogr.Geometry(osgeo.ogr.wkbLinearRing)
    # Add points
    for lon, lat in polygon:
        ring.AddPoint(lon, lat)
    # Create polygon
    poly = osgeo.ogr.Geometry(osgeo.ogr.wkbPolygon)
    poly.AssignSpatialReference(spatialReference)
    poly.AddGeometry(ring)
    # Create point
    pt = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
    pt.AssignSpatialReference(spatialReference)
    pt.SetPoint(0, point[0], point[1])
    # test
    return pt.Within(poly)

Sin embargo, no parece funcionar correctamente:

In [22]: polygon = [(-10., -10.), (-10., 10.), (10., 10.), (10., -10.), (-10., -10.)]

In [23]: point_in_polygon([0., 0.], polygon)
Out[23]: True

In [24]: point_in_polygon([359., 0.], polygon)
Out[24]: False

¿Alguna idea de lo que me falta?

5voto

Symmetric Puntos 158

Sospecho que se debe a que el subyacente GEOS La biblioteca sólo funciona en el espacio cartesiano y no en el esférico, por lo que tendrás que restar 360 a cualquier coordenada longitudinal superior a 180, lo que hace que 359 == -1. Por supuesto, todavía tendrá problemas con las características que cruzan el antimeridiano (es decir, +- 180 grados de longitud), pero puede detectar fácilmente eso y no hacer el paso de sustracción de 360.

0 votos

Gracias por su aclaración. Me llama la atención de impar que algo que tiene "geoespacial" en su nombre no cubra estas simples geoaplicaciones. ¿Alguna otra idea de cómo puedo resolver el punto-en-polígono en Python sin tener que considerar esos casos manualmente?

0voto

Developer Puntos 1191

De acuerdo con su comentario si usted está buscando el código puro en lugar de utilizar el paquete OGR, hay muchos buenos enlaces como esto: en C y esto: en Python que abordan directamente su problema por medio de códigos puros. Usted puede encontrar aún más buscando en Google;)

1 votos

Gracias, pero esos enlaces no tienen en cuenta que la Tierra es una esfera.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X