3

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.

  • Area and perimeter have nothing to do with point-in-polygon, but all three are derived from polygon coordinates. As gene says, get the coords (or geometry) of the polygon, then see if the point is inside. – Martin F Jul 27 '14 at 22:01
  • I don't have the coordinates of the polygon. All I have is the length and area. I've probably gotten myself in over my head...this is all very simple with the dataset the geodjango tutorial uses. – fromageball Jul 28 '14 at 1:33
  • If your immediate problem is "how do i get polygon coordinates from shape files" you should change it so, and omit (or downplay) anything about area, perimeter and point-in-polygon. – Martin F Jul 28 '14 at 3:55
  • Ahah! Now i see that's what you did with your new question (the one i suggested was a duplicate). Maybe delete this question... – Martin F Jul 28 '14 at 3:59
  • I thought about deleting this one, but I've gotten several pluses on it so I wasn't sure if maybe other people could use the info? – fromageball Jul 28 '14 at 11:23
8

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)
  • Thanks for the detailed response - I tried converting the format but it still doesn't work on my neighborhood data set. I was able to get the correct results using my original query with the WorldBorder data set, but I still get no results from the neighborhood data set even when I convert the format to a Point. – fromageball Jul 27 '14 at 12:05
  • Really nice answer. Is it possible to get a circular area given a point and a distance (radius) in KM? I can't find a full documentation for GEO methods. – Gocht Aug 7 '15 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.