Thursday 15 June 2017

How do I get the centroid of line with OGR


Can somebody explain how to get the centroid of a feature with ogr? I have no idea what the centroid() function itself does.


>>> linelyr = lineshp.GetLayer() 
>>> f = linelyr.GetFeature(0)

>>> g = f.geometry()
>>> g.GetPoint()
(5.85136, 49.974, 0.0) # starting point of the line

c =g.Centroid()
>>> c.GetX()
0.0

Answer



The centroid function returns the centroid of the vertices within a geometry (see the documentation). It was originally only intended for use with polygons but can now be used with other geometries. So, the start-point of a line gives no indication of where the centroid of its vertices should be and is therefore irrelevant in your example above (I'm not intending to be harsh). However, the centroid of a bunch of vertices is not the same as the mid-point of a line (i.e. halfway along it). For example, the following line string describes a square around the origin:


myWktLine = 'LINESTRING (-1 -1, -1 1, 1 1, 1 -1, -1 -1)'

myLineGeom = ogr.CreateGeometryFromWkt(myWktLine)

The mid-point traced along of the line is at (1, 1) but the centroid of this line's geometry is at (0,0):


>>> newLineGeom.Centroid().GetX()
0.0
>>> newLineGeom.Centroid().GetY()
0.0

EDIT:
To get the mid-point of the line (rather than the centroid of its vertices) you will need a function something like PostGIS' ST_LineInterpolate_Point, setting the fraction to 0.5. The same function is available in SpatiaLite if you don't want to go to the bother of installing PostGIS (and is probably available in most spatially enabled databases too).



Alternatively you could use the Shapely library, which extends the functionality of OGR nicely (without having to import your data to a database). To use Shapely, get your geometry via OGR in the way you are familiar but cast it as wkb or wkt and then use Shapely's interpolate function as follows:


import ogr
from shapely.geometry import LineString
from shapely import wkt

shapefilePath = r"C:\someFolder\someFile.shp"
ds = ogr.Open(shapefilePath)
lyr = ds.GetLayer(0)
f = lyr.GetFeature(0)
g = f.geometry().ExportToWkt()

shapelyLine = LineString(wkt.loads(g))
midPoint = shapelyLine.interpolate(shapelyLine.length/2)
print("Halfway along the line is the place where I sit:", midPoint.to_wkt())

Using the Shapely interpolate method on a linestring object constructed from 'myWktLine' in my first example correctly returns a value of 'POINT(1 1)' (when cast as wkt). Note that Shapely's interpolate function takes a distance rather than a fraction, which is why I use 'length/2' as opposed to 0.5 as you would for ST_Line_interpolate_point.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...