[GIS] Geodjango – trying to get points within a polygon

geodjangopoint-in-polygonshapefile

I'm going through the Geodjango tutorial and trying to implement in my own project – I have loaded data from a shapefile into my model/database and I have a separate model that has a PointField. I am trying to find all of the points in a polygon area but I can't figure out why my queries are returning empty.

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()

other model:

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

I have tried a few different things, including following the tutorial directions and manually typing in a Point that should definitely return something, but everything I try returns nothing.

In my view I have:

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

…and I get zero results for every WayPoint.

I have also tried(from the tutorial):

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

Nothing.

I've gone through the tutorial a few times and I can't figure out what I am doing wrong. I am new to GIS and am having trouble trying to troubleshoot the issue.

Ultimately I am trying to loop through and get all of the WayPoints in each Neighborhood, which I will eventually put on a map. But I can't even get any of my Neighborhood queries to return anything.

I feel like I must be missing something obvious here…

Edit: I went through the tutorial again and did it exactly with the WorldBorder dataset and tried my queries against that. My data set is for neighborhoods in NYC, so the point I'm using that definitely should return a polygon is a point in Times Square. When I query the neighborhood dataset I get an empty set, and I have the World Dataset at least returning back United States now. I don't understand why the neighborhood dataset is returning nothing.

Best Answer

You first need to understand the different ways to create geometries and the geometrical relations in GeoDjango (GeoDjango: GEOS API), without a database:

1) create valid geometries:

# 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) To mimic what mpoly__contains is doing, I use a list of points:

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)

Or with one list:

 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) it is the same principle with mpoly__contains=pt but you cannot mix the different geometrical formats.

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

So:

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)
Related Question