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. alt text

[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)
                if point.Within(Network.polygons[x][0].GetGeometryRef()):
    print g.GetPoint(j)

Thanks Paolo and Chris for the hints.

2 Answers 2


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.

  • grzie Paolo, but what I really need is not to compare the object references of a point with the one of a polygon, but rather the last and the first point of a linestring, which I currently access with GetPoint. There's no GetPointRef I noticed. How would I go about implementing that? Commented Jan 23, 2011 at 16:43
  • 1
    linestring is an iterable geometry collection of point (vertex), you can easily access the first and last of them.
    – capooti
    Commented Jan 23, 2011 at 22:29
  • I've tried to store the polygon objects in a list which I would later use to check against the first and last point of a linestring. However, no points are added to the list of points within a polygons for some reason: have a look here: codepad.org/Cm2BV5mp Commented Jan 25, 2011 at 11:29

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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