Measure distance between two points

[Tyler Andersen]

I'm trying to figure out how to take two Point objects and determine
the distance between them.

I have a set of Point objects, and I'd like to compute the distance
in
miles for each.

For Point objects a, b, one can "a.distance(b)", however, a number is
returned not a Distance object. I'm having trouble finding
documentation on what this value represents (I could perhaps plug it
into D(????=a.distance(b)).mi).


The docs clearly describe how work with the distance between a
reference point and points stored in a DB, but I just want the simple
operation of measuring the distance between two points, no DB
involved.

Thanks,
Tyler

[Robert Coup]

Hi Tyler,

On Mon, Jun 14, 2010 at 2:43 PM, Tyler Andersen <bigbr...@gmail.com> wrote:

I'm trying to figure out how to take two Point objects and determine
the distance between them.

I have a set of Point objects, and I'd like to compute the distance
in
miles for each.

For Point objects a, b, one can "a.distance(b)", however, a number is
returned not a Distance object. I'm having trouble finding
documentation on what this value represents (I could perhaps plug it
into D(????=a.distance(b)).mi).

In general, distances are in the units of the coordinate system of each point. So if your coordinate systems are in metres... then a.distance(b) will return the distance between a & b in metres...

>>> from django.contrib.gis import geos

>>> from django.contrib.gis.measure import D

>>>

>>> g1 = geos.Point(1000, 1000)

>>> g2 = geos.Point(2000, 1000)

>>>

>>> g1.distance(g2)

1000.0

>>> D(m=g1.distance(g2)).mi

0.62137119223733395

If you're using lat/longs for your points, you'll need to transform them into something linear (eg. an appropriate UTM zone), or do it in postgis via the GeoQuerySet.distance() method.

>>> g1 = geos.Point(176.1, -36.2, srid=4326)

>>> g2 = geos.Point(177.1, -37.2, srid=4326)

>>>

>>> g1.transform(32760) # UTM Zone 60S

>>> g2.transform(32760)

>>>

>>> g1.distance(g2)

142422.48292200203

>>> D(m=g1.distance(g2)).mi

88.497228014645742

HTH,

Rob :)

[Tyler Andersen]

On Jun 14, 1:43 pm, Robert Coup <robert.c...@koordinates.com> wrote:

Rob-

Thanks for your help. I'm a little confused here though.

If my research is correct, a UTM zone can be determined by the srid
and the lat/long (which, coincidentally are the items passed into the
constructor of a Point).




As an example, I'll take a point in SF (37.8N, 122.4W):

>>> x = Point(-122.4, 37.8, srid=4326)


...and a point further south (37.4N, 122.1W):

>>> y = Point(-122.1, 37.4, srid=4326)


Before calling transform:


>>> D(m=x.distance(y))
Distance(m=0.5)
>>> D(m=x.distance(y)).mi
0.00031068559611867048


After transform:
>>> x.transform(32610)
>>> y.transform(32610)
>>> D(m=x.distance(y)).mi
0.00031068559611867048




I guess I don't understand why transforming to an appropriate UTM zone
needs to be explicitly called out (given the lat/long/srid). On top of
that, I don't seem to be getting the answer I'm expecting. The actual
distance between x and y should be measured in tens of miles (also,
the transform seemed to have no effect).



I'm still a little confused, and I suspect I'm missing something
obvious. Any clue?


Thanks,

Tyler

[Robert Coup]

Hi Tyler,

On Tue, Jun 15, 2010 at 10:11 AM, Tyler Andersen <bigbr...@gmail.com> wrote:

If my research is correct, a UTM zone can be determined by the srid
and the lat/long (which, coincidentally are the items passed into the
constructor of a Point).

Prettymuch. For points it works nicely, assuming the points you're interestedin don't straddle a boundary. Obviously it gets trickier for lines/polygons/multipoints/etc, which is why it's not an ideal solution.

 

As an example, I'll take a point in SF (37.8N, 122.4W):

>>> x = Point(-122.4, 37.8, srid=4326)


...and a point further south (37.4N, 122.1W):

>>> y = Point(-122.1, 37.4, srid=4326)


Before calling transform:


>>> D(m=x.distance(y))
Distance(m=0.5)
>>> D(m=x.distance(y)).mi
0.00031068559611867048

So the distance is 0.5 between your two points, via Mister Pythagoras. In "degrees", which all has the usefulness of a chocolate teapot. When you create a Distance object you're saying that the number is in metres. Which it's not :( So the miles version is wrong as well.

 

After transform:
>>> x.transform(32610)
>>> y.transform(32610)
>>> D(m=x.distance(y)).mi
0.00031068559611867048

I suspect you don't have GDAL installed/working correctly. GEOSGeometry.transform() fails silently in this case (why Justin, why?!), which means you end up with what you started with (check by looking at x.wkt)

Test via:

>>> import django.contrib.gis.gdal

>>> django.contrib.gis.gdal.HAS_GDAL

True

if yours says False then transform() won't work... :(

 

Rob :)

[Tyler Andersen]

On Jun 14, 3:22 pm, Robert Coup <robert.c...@koordinates.com> wrote:

Came back, got GDAL installed. Thanks for the tip. I saw this in the
docs:

"""
While GDAL is technically not required, it is recommended. Some
features of GeoDjango (including the LayerMapping data import utility
and the geographic admin) depend on its functionality.
"""


...but I didn't think that applied here.


I see how that D() is doing The Right Thing (with respect to my
problem):


>>> from django.contrib.gis.measure import D
>>> from django.contrib.gis.geos.point import Point

>>> x = Point(-122.4, 37.8, srid=4326)

>>> y = Point(-122.1, 37.4, srid=4326)

>>> D(m=x.distance(y)).mi
0.00031068559611867048
>>> x.transform(32610)
>>> y.transform(32610)
>>> D(m=x.distance(y)).mi
32.113233025063785
>>>

[[ I do find it amusing/aggravating that transform() has different
behavior depending on whether another (unreferenced by my code)
library is installed. ]]


*however*, I'd like to bring up an earlier point:

> > If my research is correct, a UTM zone can be determined by the srid
> > and the lat/long (which, coincidentally are the items passed into the
> > constructor of a Point).
>
> Prettymuch. For points it works nicely, assuming the points you're
> interestedin don't straddle a boundary. Obviously it gets trickier for
> lines/polygons/multipoints/etc, which is why it's not an ideal solution.

I imagine there exists a function to perform this transform
automagically? I don't really follow the UTM zones to well, but it
seems to me that given a lat/long and srid, it's possible to compute
the UTM zone. It seems to me that based on this, given two Point
objects with those three data points (each), the distance should be
computable without having to add set of calls for the user (to
transform). Am I misunderstanding?



Thanks for your help,

Tyler

[Robert Coup]

On Wed, Jun 16, 2010 at 10:10 AM, Tyler Andersen <bigbr...@gmail.com> wrote:

[[ I do find it amusing/aggravating that transform() has different
behavior depending on whether another (unreferenced by my code)
library is installed. ]]

Me too. Justin? Should this be a bug? I'd certainly rather have transform explode than do Bad Things silently.

 

I imagine there exists a function to perform this transform
automagically? I don't really follow the UTM zones to well, but it
seems to me that given a lat/long and srid, it's possible to compute
the UTM zone.  It seems to me that based on this, given two Point
objects with those three data points (each), the distance should be
computable without having to add set of calls for the user (to
transform).  Am I misunderstanding?

Here's a SQL version, I doubt it's too taxing to port it to python :)

http://postgis.refractions.net/pipermail/postgis-users/2005-December/010262.html

UTM zones are 6 degrees wide, and use different coordinates within each zone. Outside the 6 degree band, they lose accuracy pretty quickly. So the distance between 

*Two* lat/long points can have their distance calculated directly without too much pain:

http://en.wikipedia.org/wiki/Great-circle_distance

Maybe something like this could be included in GeoDjango or the GEOS library...

Another solution is to use the new Geography type available in PostGIS 1.5 - it knows that "lat/longs" are on a sphere, and can calculate distances/filters/etc properly. GEOS doesn't support the geography type though so it's not available without a database...

http://www.postgis.org/documentation/manual-svn/ch04.html#PostGIS_Geography

Rob :)

[Justin Bronn]

On 6/15/10 5:52 PM, Robert Coup wrote:
>> [[ I do find it amusing/aggravating that transform() has different
>> behavior depending on whether another (unreferenced by my code)
>> library is installed. ]]
>>
>
> Me too. Justin? Should this be a bug? I'd certainly rather have transform
> explode than do Bad Things silently.

The current behavior transforms the geometry if GDAL is present and does
nothing if it isn't present. In other words, no data corruption or
other "Bad Things" occur -- the geometry remains untransformed (because
it _can't_ be, there's no GDAL available to morph the vertices).

In hindsight, an exception probably should've been raised -- but if a
ticket is opened for this 'feature', then there will be
backwards-compatibility issues to address.

-Justin

[Robert Coup]

On Thu, Jun 17, 2010 at 1:25 AM, Justin Bronn <jbr...@gmail.com> wrote:

The current behavior transforms the geometry if GDAL is present and does nothing if it isn't present.  In other words, no data corruption or other "Bad Things" occur -- the geometry remains untransformed (because it _can't_ be, there's no GDAL available to morph the vertices).

But for a method that usually mutates a geometry instance, it is bad. The caller needs to check whether the .srid attribute has changed to what you asked for... at least with transform(clone=True) you get None if it doesn't work, so code around it will explode quickly.

 

In hindsight, an exception probably should've been raised -- but if a ticket is opened for this 'feature', then there will be backwards-compatibility issues to address.

Idea 1 - issue 2 warnings?

 - a Warning that transform() is doing nothing.

 - a FutureWarning, that semantics will change from 1.5 & an error will be raised.

Idea 2 - I can't actually see any use cases where a user would want transform() to no-op unless the srids match - can anybody else? In which case I think it *could* be classed as a bug and hence ignore the deprecation policy if:

  - transform() no-ops if source & destination SRIDs match

  - we deal with geometry.srid being None to start with - i just noticed that case no-op's as well... maybe that could raise a FutureWarning too?

I guess if someone was doing the following, then they'd get different exceptions raised:

  g.transform(1234)

  if g.srid != 1234: raise MyCustomGdalAvailableError("Transform didn't work") # GDAL isn't available

  g1 = g.transform(1234, clone=True)

  if g1 is None: raise MyCustomGdalAvailableError("Transform didn't work") # GDAL isn't available

but i'd think the usual case if they were worried about it (i have this in some projects) would be:

  from django.contrib.gis.gdal import HAS_GDAL

  assert HAS_GDAL, "GDAL isn't available, we need it!"

  ... do some transform()s later ...

Ideas? Comments? Will add this to a ticket shortly.

Rob :)

[Justin Bronn]

On 6/16/10 4:17 PM, Robert Coup wrote:
> On Thu, Jun 17, 2010 at 1:25 AM, Justin Bronn<jbr...@gmail.com> wrote:
> But for a method that usually mutates a geometry instance, it is bad. The
> caller needs to check whether the .srid attribute has changed to what you
> asked for... at least with transform(clone=True) you get None if it doesn't
> work, so code around it will explode quickly.

Good points about the clone keyword and the burden placed on user for
error checking (after all frameworks should make things easier).
Someone please open a ticket for this.

Once again, an abstraction leakage is the root cause of this
side-effect. GEOS, the library, has no concept of coordinate systems
other than the integer SRID. In other words, `GEOSGeom_transform`
doesn't exist. I added `GEOSGeometry.transform` as a convenience
method, but GEOS as doesn't perform transformations, the code asks GDAL
to perform the transform -- and if GDAL isn't around then nothing
happens, which is the root cause of this thread.

(Rob knows this already, but writing it down so others may follow)

> Idea 1 - issue 2 warnings?
> - a Warning that transform() is doing nothing.

> - a FutureWarning, that semantics will change from 1.5& an error will be

> raised.
>
> Idea 2 - I can't actually see any use cases where a user would want
> transform() to no-op unless the srids match - can anybody else? In which
> case I think it *could* be classed as a bug and hence ignore the deprecation
> policy if:

> - transform() no-ops if source& destination SRIDs match

> - we deal with geometry.srid being None to start with - i just noticed
> that case no-op's as well... maybe that could raise a FutureWarning too?
>

Agreed on no-op when appropriate, and exception after if GDAL not
present. This will break backwards-compatibility, but at same time will
inform users of potential problems in their code. Trade-off seems
acceptable to me, but only if the incompatible change is thoroughly
documented. In other words, an acceptable patch will have to note this
in the `transform` docs as well as mention the backwards-incompatibility
in future 1.3 release notes.

A tougher question is whether we backport to 1.2.X.

-Justin

[Tyler Andersen]

On Jun 17, 12:24 am, Justin Bronn <jbr...@gmail.com> wrote:
> On 6/16/10 4:17 PM, Robert Coup wrote:
>
> > On Thu, Jun 17, 2010 at 1:25 AM, Justin Bronn<jbr...@gmail.com>  wrote:
> > But for a method that usually mutates a geometry instance, it is bad. The
> > caller needs to check whether the .srid attribute has changed to what you
> > asked for... at least with transform(clone=True) you get None if it doesn't
> > work, so code around it will explode quickly.
>
> Good points about the clone keyword and the burden placed on user for
> error checking (after all frameworks should make things easier).
> Someone please open a ticket for this.
>
> Once again, an abstraction leakage is the root cause of this
> side-effect.  GEOS, the library, has no concept of coordinate systems
> other than the integer SRID.  In other words, `GEOSGeom_transform`
> doesn't exist.  I added `GEOSGeometry.transform` as a convenience
> method, but GEOS as doesn't perform transformations, the code asks GDAL
> to perform the transform -- and if GDAL isn't around then nothing
> happens, which is the root cause of this thread.

Technically that was just the root cause of the end of the thread. :)

The root cause of the thread is about why the user would need to
perform additional transform() calls on two points to measure the
distance between them, when the parameter to transform() is something
of a (small, static) table lookup based on the three parameters to
point (lat, lon, srid). It seems rather inconvenient...

Thanks,

Tyler

[Robert Coup]

Hi Tyler,

On Fri, Jun 18, 2010 at 4:35 AM, Tyler Andersen <bigbr...@gmail.com> wrote:

The root cause of the thread is about why the user would need to
perform additional transform() calls on two points to measure the
distance between them, when the parameter to transform() is something
of a (small, static) table lookup based on the three parameters to
point (lat, lon, srid). It seems rather inconvenient...

We could add the haversine formula, for calculating a great-circle distance between two points. We could even find a version that uses an ellipsoid rather than a sphere so it's more accurate, but... it'd only apply for points. Is that acceptable, or just creating more weirdness & problems? Are there haversine-esque formulas for arbitrary lines/polygons/multigeometries?

I'd hope that in time the Postgis geography support code would be extricated to GEOS/elsewhere so we can transparently use it in GeoDjango.

Rob :)

[Robert Coup]

On Thu, Jun 17, 2010 at 7:24 PM, Justin Bronn <jbr...@gmail.com> wrote:

Good points about the clone keyword and the burden placed on user for error checking (after all frameworks should make things easier). Someone please open a ticket for this.

http://code.djangoproject.com/ticket/13788

Even has a patch with tests and docs.

Rob :)