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

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.

[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

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.