Point (of a linestring) within Polygon using ogr and Python

Question

I'm currently working on a project in which I need to build a topological network out of the geometry features I find in shapefiles. So far using Ben Reilly's open source project I've managed to transform linestrings into network edges, as well as detect close features (other linestrings say) and add them to the closest point so that I can run the shortest path algorithms.

But that's fine for one shapefile. However, I now need to connect features from different shapefiles into a big network graph. So for instance, if a point is within a polygon I would connect it (by connecting it I mean add a network edge - add_edge(g.GetPoint(1), g.GetPoint(2)) with the point in the next shapefile that also is within a polygon that shares a similar attribute (say, ID). Note that the polygons in the different shps only share the same IDs and not the coordinates. The points that fall within the polygons also do not share the same coordinates.

My solution to this problem was to identify the point that resides in a polygon, store it, find the point in the next shapefile that resides in the polygon with the same id and then add the networkx edge between them.

How to find if a point resides within a polygon? Well, there's a well-known algorithm: RayCasting algorithm that does that. This is where I actually got stuck though, because in order to implement the algorithm I need the coordinates of the polygon, and don't know how to access them right now, even after skimming through a documentation of the OGR's Geometry. So, the question that I'm asking is how to access the polygon points or coordinates OR is there an easier way of detecting whether a point falls within a polygon? Using python with osgeo.ogr library I coded the following:

 if g.GetGeometryType() == 3: #polygon
                c = g.GetDimension()
                x = g.GetPointCount()
                y = g.GetY()
                z = g.GetZ()

see the image for a better understanding of my issue.

[EDIT] So far I've tried storing all the polygon objects in a list with which I would then compare the linestring first and last point. But Paolo's example is related to using the Point Object reference and the Polygon Object reference, which would not work with a line object reference since not the entire line is within the polygon but rather the first or the last point of its linestring.

[EDIT3] Creating a new Geometry point object from the coordinates of the first and last point of the linestring and then using that one to compare against the polygon geometry objects saved in a list seems to work just fine:

for findex in xrange(lyr.GetFeatureCount()):
    f = lyr.GetFeature(findex)
    flddata = getfieldinfo(lyr,f,fields)
    g = f.geometry()
    if g.GetGeometryType() == 2:
        for j in xrange(g.GetPointCount()):
            if j == 0 or j == g.GetPointCount():
                point = ogr.Geometry(ogr.wkbPoint)
                point.AddPoint(g.Getx(j),g.GetY(j))
                if point.Within(Network.polygons[x][0].GetGeometryRef()):
    print g.GetPoint(j)

Thanks Paolo and Chris for the hints.

Answer 1

Shapely is cool and elegant, but why not using still ogr, with its spatial operators (in OGRGeometry class)?

sample code:

from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
polyshp = driver.Open('/home/pcorti/data/shapefile/multipoly.shp')
polylyr = polyshp.GetLayer(0)
pointshp = driver.Open('/home/pcorti/data/shapefile/point.shp')
pointlyr = pointshp.GetLayer(0)
point = pointlyr.GetNextFeature()
polygon = polylyr.GetNextFeature()
point_geom = point.GetGeometryRef()
polygon_geom = polygon.GetGeometryRef()
print point_geom.Within(polygon_geom)

Note that you need to compile GDAL with GEOS support.

Answer 2

I'm not familiar with networkx but if I understood correctly your question, you could use shapely and OGR lib for finding point in polygon from shapefile. Here is one example how it works for finding if one point (2000,1200) fails withing any polygon from one shapefile. For the result, it prints coordinates of that polygon.

from shapely.geometry import Point, Polygon
from shapely.wkb import loads
from osgeo import ogr

file1 = ogr.Open("d:\\fileWithData.shp")
layer1 = file1.GetLayerByName("fileWithData")

point1 = Point(2000,1200)

polygon1 = layer1.GetNextFeature()

while polygon1 is not None:
    geomPolygon = loads(polygon1.GetGeometryRef().ExportToWkb())
    if geomPolygon.contains(point1):
        xList,yList = geomPolygon.exterior.xy
        print xList
        print yList
    polygon1 = layer1.GetNextFeature()

I hope it helps.