Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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 networkx edges, as well as detecting close features (other linestrings say) and adding them to the closest point so that I can run shortest path algorithms.

But that's fine for one shapefile. However, I now need to connect features from different shapefiles into a big networkx graph. So for instance, if a point is within a polygon I would connect it (by connect it I mean add a networkx 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 algorith 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)
                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.

share|improve this question
add comment

2 Answers

up vote 9 down vote accepted

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.

share|improve this answer
    
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? –  bo6 Jan 23 '11 at 16:43
1  
linestring is an iterable geometry collection of point (vertex), you can easily access the first and last of them. –  capooti Jan 23 '11 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 –  bo6 Jan 25 '11 at 11:29
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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