Projection from/to the same CRS gives large Z difference since 2.3.16

[bartvanandel]

Example using EPSG:7415 (Amersfoort RD New + NAP), tested in PhantomJS

Head:

const proj4 = require('./dist/proj4.js');
proj4.defs('RD', '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +no_defs');

Using 2.3.16 or 2.3.17:

proj4('RD', 'RD', [100000, 437000, 10]);
// {
//    "0": 100000.00006519399,
//    "1": 437000.00057943614,
//    "2": -43.74852570705116
// }

Must be a regression, because 2.3.15 gives these results:

proj4('RD', 'RD', [100000, 437000, 10]);
// {
//    "0": 100000.00009487881,
//    "1": 436999.9997790383,
//    "2": 10.000132272019982
// }

Now, ideally proj4(a, a, b) should probably recognize that the "from" and "to" CRS are the same and just return b, but even without this check, I don't think such large changes in Z are acceptable. What's going on here?

[calvinmetcalf]

we really don't do anything with z coordinates, they honestly should just be passed through untouched

…

On Thu, Jan 26, 2017 at 11:25 AM Bart van Andel ***@***.***> wrote: Example using EPSG:7415 (Amersfoort RD New + NAP), tested in PhantomJS Head: const proj4 = require('./dist/proj4.js');proj4.defs('RD', '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +no_defs') Using 2.3.16 <http://releases/2.3.16>: proj4('RD', 'RD', [100000, 437000, 10]);// {// "0": 100000.00009487881,// "1": 436999.9997790383,// "2": 10.000132272019982// } Must be a regression, because 2.3.15 <http://releases/2.3.15> gives these results: proj4('RD', 'RD', [100000, 437000, 10]);// {// "0": 100000.00006519399,// "1": 437000.00057943614,// "2": -43.74852570705116// } Now, ideally proj4(a, a, b) should probably recognize that the "from" and "to" CRS are the same and just return b, but even without this check, I don't think such large changes in Z are acceptable. What's going on here? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#230>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABE4n5ymtubytC5aQXHnQ7RVESykClJlks5rWMjxgaJpZM4Lu1Ii> .

[bartvanandel]

That's essentially kind of what I expected to see, but well, that's clearly not the case here. Any ideas as to how this can happen?

[calvinmetcalf]

the datum transforms changes z values for probably a very good reason?

[bartvanandel]

Maybe, but in this example going from one point to itself (same CRS) changes the height by more than 50 meters, where previously it was more in the order of 1/10th of a micrometer. The error is about 400,000 times worse since 2.3.15 for this particular example. That can't be good.

[bartvanandel]

Just tried whether this is still the case in current master, and unfortunately, it is. Good thing I've pinned the version to 2.3.15 in our package.json file, but I'm wondering: is there any chance this issue will be addressed anytime soon? Note (as above): I'm converting the coordinate from/to the same CRS.

proj4.version
// "2.4.4-alpha"
proj4('RD', 'RD', [100000, 437000, 10]);
// {
//    "0": 100000.00006519399,
//    "1": 437000.00057943614,
//    "2": -43.74852570705116
// }

[bartvanandel]

I found the culprit. 2.3.16 introduced a refactoring where instead of modifying points in place, a new object is being created every time a modification occurs. Basically, copy-on-write. Unfortunately the z and also m properties of the input point are not always being copied.

Now, I can create a patch which copies these properties back into each newly generated point object. However, maybe it makes more sense to only apply copy-on-write for public functions and allow internal functions to modify their arguments. This results in less need for copying, meaning less checks for properties to copy, and hence, likely some performance improvements.

Anyway, to show that my (preliminary) fix is working, this is what I'm getting now, which is pretty much the same result from 2.3.15:

proj4('RD', 'RD', [100000, 437000, 10]);
// {
//    "0": 100000.00009487881,
//    "1": 436999.9997790383,
//    "2": 10.000132272019982
// }

[calvinmetcalf]

around here might be the place to do it

proj4js/lib/core.js

Line 5 in 6b0324e

function transformer(from, to, coords) {

[bartvanandel]

Currently I've changed all point = ... occurrences in lib/template.js, for instance this line only copies x and y. Isn't that the public transform function? As in aliased as proj4(from, to, coord)?

[calvinmetcalf]

no core has the public function which has the logic for allowing arrays instead of points, so we could check if the object has more then 2 keys and if it does we can copy all the other keys over at the end, I can whip up a pull.

[bartvanandel]

Okay I understand. But, as far as I can tell, there are transforms (such as the transform functions in datumUtils.js) which do convert z. So, doesn't it make more sense to make sure that the public function just creates 1 copy and have the internal functions just modify that value in place, like I suggested?

[calvinmetcalf]

ah so the reason we do the copy on write thing is because the code did a huge amount of modifications in various places making things very hard to follow, basically with any particular call you wouldn't know if it was actually modifying anything or (in a couple cases) just doing nothing so that update made the code MUCH easier to follow, as for the the datum transform that modify z, well as it turns out they don't actually work right and #192 actually complains that those values should just be passed through and that seems the least worst thing we could do.

[bartvanandel]

Can I deduce that as far as precision goes, the transforms are then only correct (up to the precision offered by the provided approximations in any formula) when z is actually zero? Coordinates with higher z could be in a different place when converted from one coordinate system to another, when projected back onto the ground (so, different x and y).

Note: we're only using Proj4 to project points onto a map. When projecting in 3D, we usually use our own transforms with full support for 3d, but that does not run in a browser.