create temporary point geometry feature from coordinates in OGR with Python

Question

I have two shapefiles containing two types of geometry features: one contains polygons and the other linestrings.

Because I'd like to use the Within method to check whether the end(or start) of a linestring falls within a polygon located in the other shapefile, I need to access the first and last points (as geometry features). At the moment I iterate through the linestring feature and can identify the points using geometry.GetPoint(x) which returns the x and y coordinates of that point.

But when I'd like to use within to check whether a given point falls in the polygon object I can't because as Paolo explained in this question I need the geometry reference (using GetGeometryRef())of a point which can be retrieved like:

...
pointshp = driver.Open('test.shp')
pointlyr = pointshp.getLayer(0)
point = pointlyr.GetNextFeature()
point_geom = point.GetGeometryRef()
point_geom.Within(polygon_geom)

but I'm actually opening a shp that has linestring features. So after skimming through Utah State University GIS slides, I came up with this:

def within(pointcoord):
driver = ogr.GetDriverByName('test.shp')
ds = driver.CreateDataSource('test.shp')
layer = ds.CreateLayer('test',geom_type=ogr.wkbPoint)
fieldDefn = ogr.FieldDefn('id',ogr.OFTInteger)
layer.CreateFiled(fieldDefn)
featureDefn = layer.GetLayerDefn()
feature = ogr.feature(featureDefn)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(pointcoord[0],pointcoord[1])
feature.SetGeometry(point)
feature.SetField('id',1)
layer.CreateFeature(feature)
pointlyr = layer.GetLayer(0)
point = pointlyr.GetNextFeature()
point_geom = point.GetGeometryRef()
for x in xrange(len(polygons))
        point_geom.Within(Network.polygons[x].GetGeometryRef())

Which looks really inefficient to me. So just for the sake of using a method(within) I need to create not only a layer, but a shapefile, when all I really need is a simple point feature from which I can get the geometry reference to compare it with the one of the polygons.

The question is, how to create a point feature (temporary one, to be destroyed) from x and y coordinates? Or is there a simpler way of doing that than the example above?

On another hand, is there any other way of getting the geometry reference of the start or end point (as a geometry feature) of a linestring feature?

Answer 1

Using the coordinates of the linestring geometry feature it is possible to simply create a geometry object like:

 point = ogr.Geometry(ogr.wkbPoint)
 point.AddPoint(g.GetX(j), g.GetY(j))#where j is the iteration in the linestring i.e. 0

that can later be used with: point.Within(polygon)

Answer 2

The Python API isn't really well documented at this time, but parallels the C/C++ APIs. 1

Saying that, you linked a non-direct source Perl API reference (at tkk.fi), may or may not be what you are looking for. However, the C++ API for the LineString (at gdal.osgeo.org/ogr/classOGRLineString.html ) has a StartPoint and EndPoint methods, so they should be reflected in Python, but they aren't because of the way that you can reference objects in Python.

So what you need to do is figure out how to create LineStrings in the python bindings, using your data. Looking at your slides, you can see that a LineString example is defined on Slide 4, but you're going to need to change that up to import it from your file.

The great thing about python, is that you can fire up the REPL, see what things are:

>>> import ogr
>>> dir(ogr.Geometry)
 ... # clipped
>>> dir(ogr.Geometry.Within)
 ... # clipped
>>> print ogr.Geometry.Within
<unbound method Geometry.Within>
>>> print ogr.Geometry.Within.__doc__
Within(self, Geometry other) -> bool

int OGR_G_Within(OGRGeometryH hThis,
OGRGeometryH hOther) 
>>> line = ogr.Geometry(ogr.wkbLineString) # the constant value in the constructor here determines the type of geometry it is
>>> print line
LINESTRING EMPTY
>>> line.AddPoint(10,10)
>>> line.AddPoint(20,20)
>>> line.AddPoint(30,30)
>>> line.GetPoint(1)
(20.0, 20.0, 0.0)
>>> line.GetPoint(1) # start point
(20.0, 20.0, 0.0)
>>> line.GetPoint(line.GetPointCount()-1) # end point
(30.0, 30.0, 0.0)

Extrapolating a bit, you can just call Within on the start and end points on the line for each polygon object.