Friday, 3 March 2017

Get LineString length in meters (Python, GEODjango)


In a Django model, I have a LineString field which holds GPS coordinates of a street route. I know that I the .length holds the length of the LineString, but I don't know what unit the returned value is in.


The only docs I've found (https://docs.djangoproject.com/en/1.9/ref/contrib/gis/geos/#django.contrib.gis.geos.GEOSGeometry.length) don't say that.


I tried with simple lines and got the following results:


>>> l = LineString([(0, 0), (0, 1), (1, 1)])
>>> l.length
2.0


>>> l = LineString([(0, 0), (360, 0)])
>>> l.length
360.0

the second one, in my app, should be Earth's equatorial circumference expressed in meters (instead of 360).


So I guess the length is "unit-less", which leads to my question:



  1. how can I convert the .length to meters, if the coordinates are GPS coordinates, expressed in degrees?

  2. (less important, I'm just curious) the unit of measure of the returned value.




Answer



The length attribute of a LineString will return a value in the units of the coordinate system of the geometry. The coordinate system in your case is not specified, so there is not much meaning to the length.


Your input of 360 is in degrees, so you assume an unprojected coordinate system, such as WGS84. This is why you get a length in degrees, even without specifying a coordinate system.


In order to get a length in meters, you have two options:



  1. Specify a coordinate system on the geometry and then transform it into a projected coordinate system.

  2. Input your coordinates with already projected values that have the unit of meters, and specify the corresponding coordinate system.


Assuming you have your data in unprojected (lat/lon) coordinates, the first option is probably easier to implement. Using your example:


# Create a line string

>>> line = LineString([(0, 0), (360, 0)])

# Specify srid using [WGS84][1]
>>> line.srid = 4326

# Transform into projected coordinate system (using web mercator)
>>> line.transform(3857)

# Line length. This is zero within numerical precision, because the start
# point (0, 0) is the same as the end point (0, 360)) on a map.

>>> line.length
>>> 1.1329847282581795e-08

# An example forcing the line through the other side of the world,
# and specifying srid on definition.
>>> line = LineString([(0, 0), (180, 0), (360, 0)], srid=4326)
>>> line.tansform(3857)
>>> line.length
>>> 40075016.6855784 # Close to the size of the equator


# To illustrate the dependency on the coordinate system, use an
# European centred coordinate system (LAEA).
>>> line.transform(3035)
>>> line.length
>>> 34119703.74667985

Links to coordinate system definitions: WGS84, Pseudo Mercator (aka Web Mercator), LAEA.


A final sidenote: the web mercator projection is actually not a very good example, as it mixes both concepts of projected vs unprojected, but its used often in web mapping and gives a number close to the size of the equator in your example. See https://en.wikipedia.org/wiki/Web_Mercator


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