Geodjango - tratando de obtener puntos dentro de un polígono

Question

Estoy repasando el tutorial de Geodjango y tratando de implementarlo en mi propio proyecto - He cargado datos de un shapefile en mi modelo/base de datos y tengo un modelo separado que tiene un PointField. Estoy tratando de encontrar todos los puntos en un área de polígono, pero no puedo entender por qué mis consultas están regresando vacío.

class Neighborhood(models.Model):
    boroname = models.CharField(max_length=13)
    borocode = models.IntegerField()
    countyFIPS = models.CharField(max_length=3)
    nta_code = models.CharField(max_length=4)
    nta_name = models.CharField(max_length=75)
    shape_leng = models.FloatField()
    shape_area = models.FloatField()
    mpoly = models.MultiPolygonField(srid=4326)
    objects = models.GeoManager()   

otro modelo:

class WayPoint(models.Model):
    complaint = models.ForeignKey('Complaint')
    geometry = models.PointField(srid=4326)
    objects = models.GeoManager()

He intentado algunas cosas diferentes, incluso siguiendo las instrucciones del tutorial y escribiendo manualmente un punto que definitivamente debería devolver algo, pero todo lo que intento no devuelve nada.

En mi opinión sí:

points = WayPoint.objects.all()
for point in points:
    pt = point.geometry #so this gets the PointField
    pts = Neighborhood.objects.filter(mpoly__contains=pt)

...y obtengo cero resultados para cada WayPoint.

También he probado (desde el tutorial):

pnt_wkt = 'POINT(-73.3385 40.7245)'
obj = Neighborhood.objects.filter(mpoly__contains=pnt_wkt)

Nada.

He repasado el tutorial varias veces y no consigo saber qué estoy haciendo mal. Soy nuevo en el SIG y estoy teniendo problemas para tratar de resolver el problema.

En última instancia, estoy tratando de hacer un bucle y obtener todos los WayPoints en cada Barrio, que eventualmente pondré en un mapa. Pero ni siquiera puedo conseguir que ninguna de mis consultas de barrio devuelva nada.

Siento que debo estar perdiendo algo obvio aquí...

Edición: He vuelto a leer el tutorial y lo he hecho exactamente con el conjunto de datos WorldBorder y he probado mis consultas con él. Mi conjunto de datos es para los barrios de Nueva York, por lo que el punto que estoy usando que definitivamente debe devolver un polígono es un punto en Times Square. Cuando consulto el conjunto de datos de los barrios, obtengo un conjunto vacío, y el conjunto de datos mundial me devuelve al menos Estados Unidos. No entiendo por qué el conjunto de datos del barrio no devuelve nada.

Answer 1

Primero hay que entender las diferentes formas de crear geometrías y las relaciones geométricas en GeoDjango ( GeoDjango: API GEOS ), sin base de datos:

1) crear geometrías válidas:

# with Point, Polygon objects of GeoDjango
from django.contrib.gis.geos import Point, Polygon, 
poly = Polygon(((0.0, 0.0), (0.0, 50.0), (50.0, 50.0), (50.0, 0.0), (0.0, 0.0)))
point = Point(10,10)
print poly.contains(point)
True
# with wkt string format
from django.contrib.gis.geos import fromstr
poly  = fromstr('POLYGON ((0 0, 0 50, 50 50, 50 0, 0 0))')
point1 = fromstr('Point(10,10)')
print poly.contains(point1)
True
# with GEOSGeometry objects (wkt and with srs here)
from django.contrib.gis.geos import GEOSGeometry
poly = GEOSGeometry('SRID=32140;POLYGON ((0 0, 0 50, 50 50, 50 0, 0 0))')
point1 = GEOSGeometry('SRID=32140;POINT(10 10)')
print poly.contains(point1)
True
# and 
point2 = GEOSGeometry('SRID=32140;POINT(50 50)')
print poly.contains(point2)
False

2) Para imitar lo que mpoly__contains está haciendo, uso una lista de puntos:

points=[points1, points2,GEOSGeometry('SRID=32140;POINT(35 45)'),GEOSGeometry('SRID=32140;POINT(55 45)')]
for point in points:
    if poly.contains(point):
       print point
POINT (10.00 10.00)
POINT (35.00 45.00)

O con una lista:

 points_contained = [point.json for point in points if poly.contains(point)]
 print points_contained
 [u'{ "type": "Point", "coordinates": [ 10.0, 10.0 ] }', u'{ "type": "Point", "coordinates": [ 35.0, 45.0 ] }']

3) es el mismo principio con mpoly__contains=pt pero no se pueden mezclar los diferentes formatos geométricos.

pnt_wkt = 'POINT(-73.3385 40.7245)' 
print type(pnt_wkt)
<type 'str'> #-> a simple string
pnt = fromstr(pnt_wkt) 
type(pnt)
<class 'django.contrib.gis.geos.point.Point' # -> a GeoDjango Point

Así que:

from django.contrib.gis.geos import Point
point = Point(-73.3385 40.7245) 
# or 
point = fromstr('POINT(-73.3385 40.7245)')
# and
obj  = WorldBorder.objects.get(mpoly__intersects=point)