openlayers / openlayers

OpenLayers
https://openlayers.org
BSD 2-Clause "Simplified" License
11.45k stars 3.04k forks source link

linestring getLength not returning expected value #3533

Closed alexandre-melard closed 5 years ago

alexandre-melard commented 9 years ago

I am using getLength to retrieve the linestring length.

For the same segment: 1- when using google map measure tool, I get 228m 2- when using IGN geoportail measure tool, I get 228m 3- when I use e.feature.getGeometry().getLength() I get 330m

Here are the flat coordinates: e.feature.getGeometry().getFlatCoordinates() : [571382.4214041593, 5723486.068714521, 571593.8175605105, 5723741.65502785]

in 4326: [5.132815622245775, 45.644023326845485, 5.134714626228319, 45.64562844964627]

When I check the coordinates position on either ol3 or google map, I get the same points. The difference must come from the calcul...

Did I miss something and I should not use the getLength method? Please give me some direction if you think this is not an issue.

alexandre-melard commented 9 years ago

I've done many tests and I get a 1.44 factor all the time between the distances. Maybe it's a unit problem. But I though the default unit in ol3 was meter...

alexandre-melard commented 9 years ago

well for the moment, I'll just divide my result by 1.44

bartvde commented 9 years ago

I am not sure but maybe the other frameworks use geodesic measurements? Check the measure example in ol3 which has a geodesic option.

bartvde commented 9 years ago

Please don't cross-post in multiple systems:

http://gis.stackexchange.com/questions/142062/openlayers-3-linestring-getlength-not-returning-expected-value

bill-chadwick commented 9 years ago

I think Your coordinate spacing is too small for geodesic to make a difference. Do you get a better result with a latitude near the equator?

alexandre-melard commented 9 years ago

@bartvde Crossposting is a good way to rush people, I don't expect everybody to read the ol3 issues.

I've changed my length calcul to use geodesic, I get an exact answer now.

    var wgs84Sphere = new ol.Sphere(6378137);
    var coordinates = e.feature.getGeometry().getCoordinates();
    var length = 0;
    var sourceProj = map.getView().getProjection();
    for (var i = 0, ii = coordinates.length - 1; i < ii; ++i) {
      var c1 = ol.proj.transform(coordinates[i], sourceProj, 'EPSG:4326');
      var c2 = ol.proj.transform(coordinates[i + 1], sourceProj, 'EPSG:4326');
      length += wgs84Sphere.haversineDistance(c1, c2);
    }

Now, I don't understand why the getLength doesn't return the same value. And I am also a bit worried by the fact that all ol.Sphere functions (even the constructor) are experimental. (see http://openlayers.org/en/v3.4.0/apidoc/ol.Sphere.html?unstable=true)

bartvde commented 9 years ago

No cross-posting is actually quite rude, and increases your chances of getting ignored. Most developers use github notifications and are subscribed to the openlayers-3 tag on stack overflow, they will get notified twice.

Please only open up a github issue once you've confirmed and are sure it's a bug in the library.

alexandre-melard commented 9 years ago

Ok, fine, I was not aware of your organisation, so is it better to post on stackexchange first, and then post an issue here if the bug is confirmed? I didn't mean to be rude but just to help by alerting on a suspect behavior

marcjansen commented 9 years ago

@mylen Yes, questions go to stackoverflow, and if a bug / misbehaviour is confirmed there, open an issue here. If you are quite sure you have discovered a bug, you may skip the stackoverflow question, but of course it wouldn't hurt to ask there first.

bill-chadwick commented 9 years ago

Mylen, you get the factor of 1.4 on EPSG:3857 (Web Mercator) because 1 / cos (latitude of 45.6 N) is about that and your points are around 45.6 degrees North.

With Web Mercator, scale changes with latitude (see the Wikipedia article on the Mercator Projection for more details). The method ol.geom.LineString.prototype.getLength will only work (give accurate answers in metres) for short lines near the equator.

Currently, the comment / API doc for the method is this

/**

I think that is not telling the whole story - IMHO a stronger health warning is needed. It should perhaps say that the returned units are not always metres. There are other object methods in ol.geom that will have similar issues e.g. Polygon.getArea

Perhaps this issue should be reopened until the docs are sorted.

Personally I am not convinced that the methods should be 'api stable' as they are so dangerous for use with Web Merctor. Even relative comparisons of lengths / areas are not possible with geometries at significantly different latitudes.

marcjansen commented 9 years ago

Pull requests are always welcome.

bill-chadwick commented 9 years ago

I could do that. Its a fair fit with my pending graticule work. On 10 Apr 2015 16:42, "Marc Jansen" notifications@github.com wrote:

Pull requests are always welcome.

— Reply to this email directly or view it on GitHub https://github.com/openlayers/ol3/issues/3533#issuecomment-91596900.

probins commented 9 years ago

I am not convinced that the methods should be 'api stable' as they are so dangerous for use with Web Merctor

I think the assumption is (and IMO should be) that those using the API know the limitations of using global projections like Mercator. The method is 'stable' in that it's not going to change. Whether it's useful or not is another matter :-) IMO OL2 was more helpful in this respect in having getGeodesicLength/Area methods (Vincenty rather than Haversine) as well as simple planar ones. Haversine is exported currently. Vincenty is available in the source (under ellipsoid/), but this is not currently exported in the standard build.

alexandre-melard commented 9 years ago

@probins, when you make an assumption :

the assumption is (and IMO should be) that those using the API know the limitations of using global projections like Mercator

You mean that the usage of openlayers is targeted for expert? I'm not sure that only geographers are using it . If there is a problem with length calculation in the Mercator projection, why is it used as default on the map, could you give me an example of a use of getLength in a mercator projection?

bill-chadwick commented 9 years ago

I agree with Mylen. More of a worry to me is that there is no good/easy solution.

The Geometry classes do not record the projection of their points, so that length and area methods returning answers in metres and square metres (or other projection units) can not be implemented by them without significant changes. Even if they did know their own projection, we would not want per projection calculations in the geometry classes. IMHO the length / area maths belongs in the projection classes.

Either the length and area methods in ol.geom.Geometry classes would have to take an ol.proj.ProjectionLike parameter so that appropriate maths can be done, or else the area and length methods are moved elsewhere - e.g. to the ol.proj.Projection objects.

The latter make most sense to me. An ol.proj. Projection could be required to provide functions to calculate a Geometry length / perimeter, in projection units (often metres), and its area, in square projection units (often square metres). Having to pass a ProjectionLike into ol.geom.Geometry class length and area methods when the developer already knows (or should) the projection of the coordinates in the Geometry, would seem odd.

For 'local' projections like UTM and the UK National Grid (which have very near uniform scale across their extent) , simple maths as already exists in ol.geom.flat.lengthflatgeom and ol.geom.flat.areaflatgeom, will suffice. Perhaps some of the old OL2 code could be reused for the 'global' Web Mercator (3857) and Lat/Lon (4326) length and area calculations.

probins commented 9 years ago

but won't the methods in ol.Sphere/ol.Ellipsoid, as used in the measure example, suffice? This is independent of projection. Do you actually need per projection calculations? I'd agree that it's potentially confusing for people who aren't aware of these issues (99+% of the population) that getLength isn't appropriate for 3857, but I don't think it's the job of API docs to discuss these issues. They could usefully point people to the measure example when using 3857 and similar projections. And, as I say, IMO the API would be better if there were (stable) geodesic methods on the geometries as in OL2.

ahocevar commented 9 years ago

Re-opening because this is turning into a good discussion.

bill-chadwick commented 9 years ago

Thanks for re-opening. There are several solution options to think about.

probins commented 9 years ago

the text in the measure example could be improved too. "Earth curvature is not taken into account" isn't really correct. The whole point of projections is that they do take curvature into account. The issue is that they have different priorities. Mercator isn't intended for measuring distance or area, so getLength/Area isn't appropriate here. On the other hand, I would be pretty sure that getLength, for example, on the OSGB projection that @bill-chadwick mentions will be more accurate than Haversine.

bill-chadwick commented 9 years ago

I am coming round to the idea that the Geometry classes should know the projection of their coordinates. We could add an optional projectionLike parameter to the Geometry constructors (defaulting to EPSG:3857) and store the projectionLike as a class member. This new member would be updated if the Geometry is transformed using ol.geom.Geometry.prototype.transform but would need a setter method call after using ol.geom.Geometry.prototype.applyTransform (or the applyTransform method should be removed).

If the Geometry classes know their projection, then the getLength and getArea methods could use ol.Sphere / ol.Ellipsoid methods for projections where ol.proj.Projection.prototype.isGlobal returns true and the existing ol.geom.flat.lengthflatgeom and ol.geom.flat.areaflatgeom if isGlobal returns false. The returned length units would normally be metres and the area units metres x metres. However for non-global projections with other distance units, like feet, then the returns would be in e.g. feet or feet x feet.

ahocevar commented 9 years ago

We have discussed to let geometries know their projection in the past, but decided against it for several reasons.

probins commented 9 years ago

alternatively, you could move that decision into userland, and pass in a parameter: say, projected plane as default, or 'h' for Haversine, 'v' for Vincenty: getLength('h')

Another oddity of getLength at the moment is that if your view is in 4326 you get a 'length' in decimal degrees - I struggle to think how you could use that :-)

bill-chadwick commented 9 years ago

Vincenty and Haversine work on 4326 coords, how would we know the projection to transform from ?

probins commented 9 years ago

true. 'Pass in 2 parameters ...'

ahocevar commented 9 years ago

true. 'Pass in 2 parameters ...'

Makes sense to me, and would be consistent with other functions.

bill-chadwick commented 9 years ago

Should the params should be optional and the defaults be Haversine and 3857 ?

bill-chadwick commented 9 years ago

or perhaps default to Haversine for 3857 and 4326 and Projected Plane for non 'global' projections.

ahocevar commented 9 years ago

Haversine as default yes, but no default for the projection. So projection mandatory, and method optional.

bill-chadwick commented 9 years ago

I did a bit on this last night and plan to have a PR ready in a few days. I will add some new length and area method tests.

I propose calculation methods of 'CARTESIAN', 'HAVERSINE' and 'VINCENTY' though 'CARTESIAN', 'SPHERICAL' and 'ELLIPSOIDAL' might do instead (any opinions as to which are best?). The functions do not need a projection for the 'CARTESIAN' method though they could return NaN if projection units are degrees (4326) and the method is 'CARTESIAN'.

Calculation method will default to 'HAVERSINE' for global projections and 'CARTESIAN' for non global projections.

probins commented 9 years ago

another issue here is that of units. A really friendly function would give you the ability to define the units 'give me the length in miles'. That could be kept in userland though. Units can also be 'pixels' in which case only cartesian.

IMO 'haversine' and 'vincenty' is better than spherical/ellipsoid, as it's more specific. And I would prefer an abbreviated form like 'h' and 'v' - because I don't like typing :-)

probins commented 9 years ago

no default for the projection. So projection mandatory

this will always be the view projection though, won't it? Perhaps add a note in the docs that this will normally be view.getProjection()

bill-chadwick commented 9 years ago

I guess we would want an enum of geometry distance units and a default units parameter value of metres. We can't really reuse the ol.proj.Units enum as it includes degrees. We would have acres, hectares etc for geometry area units.

bill-chadwick commented 9 years ago

To avoid bloat, I think perhaps the length and area methods should always return metres and metres x metres, with any unit conversions being done by the user as required. A unit parameter could always be added later if deemed necessary.

probins commented 9 years ago

projection mandatory

that also makes this a breaking change, so there'll have to be release notes for this too.

probins commented 9 years ago

the text in the measure example could be improved

I was going to submit a PR for this, but the example should be changed to reflect this new method anyway (it will be simpler), so I'll hold off on this.

bill-chadwick commented 9 years ago

AARGHH!

For convenience, I would like to build ol3 on Windows but the build system is so fragile it drives me mad.

There is some issue with generate-info.js and ('child_process').spawn that is very unreliable. Sometimes (occasionally) the build works other (most) times it fails without generating info.json. A fix would be MOST welcome.

Incidentally we will need to look at / update ol.interaction.Draw as it calls LineString.getLength

ahocevar commented 9 years ago

Improvements to the build system are an ongoing effort and should land in master soon. As a result, building on Windows should be possible again soon.

bill-chadwick commented 9 years ago

Looks like ol.interaction.Draw is going to unravel my plans a bit.

ol.interaction.Draw uses geometry.setRadius(sketchLineGeom.getLength()); to set the radius of a circle. We can use event.map.getView().getProjection() to find the projection to use for an updated lineString.getLength method in ol.interaction.Draw.prototype.modifyDrawing_ but ...

For the draw interaction to work on a 4326 View I think we would have to allow the use of the Cartesian length calculation method for the 4326 projection (and other projections with units of degrees). To draw circles with non degrees projections requires that lineString.getLength must always return its result in projection units.

bill-chadwick commented 9 years ago

http://www.andre-simon.de/doku/ansifilter/en/ansifilter.php is useful for looking at test output on Windows.

bill-chadwick commented 9 years ago

I have committed the getLength changes to a branch here

https://github.com/bill-chadwick/ol3/commit/49ca11198a0a627e06722c0bd648761f1b88184e

the getArea change is still TODO.

I have added several tests for the updated getLength method

Feedback welcome.

To think about some more ...

4326 geometry simplification (Doouglas Peuker) at extreme latitudes (where x coords are smaller than y coords)

ol.geom.Circle.setRadius - 'The radius is in the units of the projection.' is a misleading comment for 3857

bill-chadwick commented 9 years ago

For 3857, we can calculate LineString length as a summation of line segment rhumbline distances (on the 3857 sphere). I believe this should be the default method for 3857 as it is the length of the lines as plotted on the map, rather than a summation of great circle segment length.

see

http://en.wikipedia.org/wiki/Mercator_projection - Formulae for distance section.

This should have better performance than using Haversine or Vincenty.

I think the maths for rhumb line length should go in epsg3857projection.js

probins commented 9 years ago

I'm not sure this is worth it: I think these days build size and latency are more the issue than performing calculations. And an advantage of having Haversine as default is it should produce the same result as measuring the same line on e.g. GMaps.

bill-chadwick commented 9 years ago

I will document that Haversine and Vincenty are great circle methods. Calculation of Rhumb line (loxodrome) length from lat/lon and EPSG:3857 coords uses very similar maths. So for completeness, I do plan to add another Geometry Length Calculation method for Rhumb Lines on the sphere, though will leave out Ellipsoidal Rhumb line maths for now.

see http://www.movable-type.co.uk/scripts/latlong.html

which I plan to use for data for the tests.

The existing measure example at

http://openlayers.org/en/v3.4.0/examples/measure.html?q=measure

is really rather wrong in several respects.

bill-chadwick commented 9 years ago

I plan to add a rhumb line length method to ol.Sphere where there is currently the comment

// FIXME add rhumb lines

bill-chadwick commented 9 years ago

I have added spherical Rhumb Line as a length calculation method. I now need to add rhumb line tests for sphere.js where the new rhumb length maths is.

bill-chadwick commented 9 years ago

PR here https://github.com/openlayers/ol3/pull/3609

bill-chadwick commented 9 years ago

Updated PR here https://github.com/openlayers/ol3/pull/3609

stale[bot] commented 5 years ago

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.