New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] 4D API cs2cs style #731

Merged
merged 4 commits into from Feb 1, 2018

Conversation

Projects
2 participants
@busstoptaktik
Member

busstoptaktik commented Jan 6, 2018

This PR is in response to issue #700. It is WIP and currently untested, but I push the PR in order to check (or actually re-establish) orthogonality to pj_transform and cs2cs.

UPDATE: Now very much tested and very little WIP. Getting ready for merging.

@busstoptaktik

This comment has been minimized.

Show comment
Hide comment
@busstoptaktik

busstoptaktik Jan 6, 2018

Member

Neat - we now have (to the extent test coverage allows) confirmed the orthogonality of this PR to the cs2cs and pj_transform sub-systems.

Next step will be to construct test data, to confirm (or probably rather reject) the hypothesis that the code does what it's supposed to. Despite just having to follow the recipe by @kbevers from #700, I doubt that I have managed to get everyting right at first shot.

Member

busstoptaktik commented Jan 6, 2018

Neat - we now have (to the extent test coverage allows) confirmed the orthogonality of this PR to the cs2cs and pj_transform sub-systems.

Next step will be to construct test data, to confirm (or probably rather reject) the hypothesis that the code does what it's supposed to. Despite just having to follow the recipe by @kbevers from #700, I doubt that I have managed to get everyting right at first shot.

@kbevers

Just a few comments. Seems like it's going in the right direction.

Show outdated Hide outdated src/pj_fwd.c Outdated
Show outdated Hide outdated src/PJ_helmert.c Outdated
Show outdated Hide outdated src/pj_fwd.c Outdated
Show outdated Hide outdated src/pj_fwd.c Outdated
Show outdated Hide outdated src/pj_fwd.c Outdated
Show outdated Hide outdated src/proj_4D_api.c Outdated
These tests are mostly based on the same material as those in
more_builtins.gie, since we are testing the same kinds of things,
but provided through a different interface.

This comment has been minimized.

@kbevers

kbevers Jan 9, 2018

Member

This should also contain a bunch of tests using the epsg init-file. Test data is easy to prepare with cs2cs. I can come up with some suggestions if needed.

@kbevers

kbevers Jan 9, 2018

Member

This should also contain a bunch of tests using the epsg init-file. Test data is easy to prepare with cs2cs. I can come up with some suggestions if needed.

This comment has been minimized.

@busstoptaktik

busstoptaktik Jan 9, 2018

Member

That would be nice. What is really needed (and will probably cause hiccups) is cases where more than one of the cs2cs style options are used.

@busstoptaktik

busstoptaktik Jan 9, 2018

Member

That would be nice. What is really needed (and will probably cause hiccups) is cases where more than one of the cs2cs style options are used.

This comment has been minimized.

@kbevers

kbevers Jan 9, 2018

Member

I'll merge #705 s soon as all the tests pass. Sync your branch with master and you'll have a bunch of test material in test/gigs/. I think that should be enough for this purpose.

@kbevers

kbevers Jan 9, 2018

Member

I'll merge #705 s soon as all the tests pass. Sync your branch with master and you'll have a bunch of test material in test/gigs/. I think that should be enough for this purpose.

This comment has been minimized.

@busstoptaktik

busstoptaktik Jan 13, 2018

Member

Have sync'ed, but currently drowning in all the "almost OK" tests. I think there is a bug in the interpretation of the GIGS tolerance parameter.

If I understand @micahcochran's comments correctly, they are considered to be max absolute coordinatewise deviations (i.e. in decimal degrees for latitude/longitude pairs), whereas in the .gie version, they are translated into "same numerical value but interpreted as having a meter unit".

This makes a.o. all NADCON tests fail, because 1e-7 m (0.0001 mm) is much less than 1e-7 degree (1 cm at the Equator), at least if you are more than a few minutes of arc away from the poles.

@busstoptaktik

busstoptaktik Jan 13, 2018

Member

Have sync'ed, but currently drowning in all the "almost OK" tests. I think there is a bug in the interpretation of the GIGS tolerance parameter.

If I understand @micahcochran's comments correctly, they are considered to be max absolute coordinatewise deviations (i.e. in decimal degrees for latitude/longitude pairs), whereas in the .gie version, they are translated into "same numerical value but interpreted as having a meter unit".

This makes a.o. all NADCON tests fail, because 1e-7 m (0.0001 mm) is much less than 1e-7 degree (1 cm at the Equator), at least if you are more than a few minutes of arc away from the poles.

This comment has been minimized.

@kbevers

kbevers Jan 13, 2018

Member

Which tests fail? There's a few of them that doesn't work in the current python setup as well. It might be that it is only those tests that fail, in which case you can disregard them, since the problems are not related to the work you are doing here.

@kbevers

kbevers Jan 13, 2018

Member

Which tests fail? There's a few of them that doesn't work in the current python setup as well. It might be that it is only those tests that fail, in which case you can disregard them, since the problems are not related to the work you are doing here.

This comment has been minimized.

@kbevers

kbevers Jan 26, 2018

Member

This makes pj_fwd/pj_inv work the same way as pj_transform: omitting dual identical helmert shifts, so it doesn't really change the test - it rewrites the test in order to act identically to the json case.

I understand why you've done it, I just don't think it is going to work. The end goal of this PR is, amongst other things, to be able to put to epsg CRS definitions into proj_create_crs_to_crs() and get correct results when using the resulting PJ in transformations. While your solution fixes the test, it doesn't fix the transformation in the test when used in the wild (in the 4D API).

@kbevers

kbevers Jan 26, 2018

Member

This makes pj_fwd/pj_inv work the same way as pj_transform: omitting dual identical helmert shifts, so it doesn't really change the test - it rewrites the test in order to act identically to the json case.

I understand why you've done it, I just don't think it is going to work. The end goal of this PR is, amongst other things, to be able to put to epsg CRS definitions into proj_create_crs_to_crs() and get correct results when using the resulting PJ in transformations. While your solution fixes the test, it doesn't fix the transformation in the test when used in the wild (in the 4D API).

This comment has been minimized.

@kbevers

kbevers Jan 26, 2018

Member

In that context, please note that I have merged your previous PR over at busstoptaktik/PROJ.4.

Duly noted.

@kbevers

kbevers Jan 26, 2018

Member

In that context, please note that I have merged your previous PR over at busstoptaktik/PROJ.4.

Duly noted.

This comment has been minimized.

@busstoptaktik

busstoptaktik Jan 26, 2018

Member

The end goal of this PR is, amongst other things, to be able to put to epsg CRS definitions into proj_create_crs_to_crs() and get correct results when using the resulting PJ in transformations

proj_create_crs_to_crs() is a much different situation:

First, we do not necessarily need this, if we can live with the time and slight precision penalty: For typical crs_to_crs uses we do just one half-roundtrip, and do not accumulate helmert imprecision.

Second, we can implement the pj_transform trick as an optimization in a much simpler way than in the generic pipeline case: Here we know we have a one input, one output situation, and if both are from epsg definitions we can check whether the towgs84=... incantations are identical and bypass Helmert by setting the towgs84=0.0.0 thingy at the proj=pipeline level.

This can be economically implemented by first changing

    strcpy(buffer, "+proj=pipeline +step +init=");

to

strcpy(buffer, "+proj=pipeline +towgs84=0,0,0 +step +init=");

Then going on with the setup, but adding some 15 lines of code, checking the case of identical Helmert parameters (instantiate, check), and then finally:

if (helmerts_not_identical)
    return P;
proj_destroy (P);
t = strstr (buffer, "towgs84");
strcpy (t, "tonight"); /* or xxxxxxx or whatever fits */
P = proj_create(ctx, buffer);
return P;

Simpler and more localized than solving the general problem in pj_pipeline

So, by and large, see the handheld introduction of towgs84=0,0,0 in the *.gie files as a temporary exercise, helping pinpoint some cases where we can do somewhat mechanically introduced optimizations in a comming (version 6-ish) early version of a late binding EPSG implementation.

@busstoptaktik

busstoptaktik Jan 26, 2018

Member

The end goal of this PR is, amongst other things, to be able to put to epsg CRS definitions into proj_create_crs_to_crs() and get correct results when using the resulting PJ in transformations

proj_create_crs_to_crs() is a much different situation:

First, we do not necessarily need this, if we can live with the time and slight precision penalty: For typical crs_to_crs uses we do just one half-roundtrip, and do not accumulate helmert imprecision.

Second, we can implement the pj_transform trick as an optimization in a much simpler way than in the generic pipeline case: Here we know we have a one input, one output situation, and if both are from epsg definitions we can check whether the towgs84=... incantations are identical and bypass Helmert by setting the towgs84=0.0.0 thingy at the proj=pipeline level.

This can be economically implemented by first changing

    strcpy(buffer, "+proj=pipeline +step +init=");

to

strcpy(buffer, "+proj=pipeline +towgs84=0,0,0 +step +init=");

Then going on with the setup, but adding some 15 lines of code, checking the case of identical Helmert parameters (instantiate, check), and then finally:

if (helmerts_not_identical)
    return P;
proj_destroy (P);
t = strstr (buffer, "towgs84");
strcpy (t, "tonight"); /* or xxxxxxx or whatever fits */
P = proj_create(ctx, buffer);
return P;

Simpler and more localized than solving the general problem in pj_pipeline

So, by and large, see the handheld introduction of towgs84=0,0,0 in the *.gie files as a temporary exercise, helping pinpoint some cases where we can do somewhat mechanically introduced optimizations in a comming (version 6-ish) early version of a late binding EPSG implementation.

This comment has been minimized.

@kbevers

kbevers Jan 26, 2018

Member

Okay, that works for me. Do you want to make that change to proj_create_crs_to_crs? It fits nicely within this PR.

@kbevers

kbevers Jan 26, 2018

Member

Okay, that works for me. Do you want to make that change to proj_create_crs_to_crs? It fits nicely within this PR.

This comment has been minimized.

@busstoptaktik

busstoptaktik Jan 26, 2018

Member

Sure, but I'm offline until Monday, so not until then.

@busstoptaktik

busstoptaktik Jan 26, 2018

Member

Sure, but I'm offline until Monday, so not until then.

Show outdated Hide outdated src/pj_fwd.c Outdated

@busstoptaktik busstoptaktik self-assigned this Jan 13, 2018

@busstoptaktik busstoptaktik modified the milestones: 4.9.4, 4.10.0 Jan 13, 2018

@kbevers

This comment has been minimized.

Show comment
Hide comment
@kbevers

kbevers Jan 26, 2018

Member

Regarding test 5208 (and possibly others). The problem is due to prime meridian not being adjusted correctly in the inverse operation, as examplified here:

λ echo 0 0 | bin\cct.exe +proj=latlong +pm=paris
 -2.3372291667    0.0000000000        1.#INF        1.#INF

λ echo 0 0 | bin\cct.exe -I +proj=latlong +pm=paris
  0.0000000000    0.0000000000        1.#INF        1.#INF # this should be  2.3372291667    0.0000000000

The forward operation is handled correctly because this

    case PJ_IO_UNITS_RADIANS:
        if (INPUT_UNITS==PJ_IO_UNITS_RADIANS)
            break;

happens in pj_fwd_finalize (line 131 and below), which basically omits the adjusted the longitude back to the greenwhich meridian. A similar construct is missing in pj_inv_prepare (remember we are going backwards now), so the longigude is adjusted to the desired prime meridan and then back to greenwich. The last bit breaks the test.

Member

kbevers commented Jan 26, 2018

Regarding test 5208 (and possibly others). The problem is due to prime meridian not being adjusted correctly in the inverse operation, as examplified here:

λ echo 0 0 | bin\cct.exe +proj=latlong +pm=paris
 -2.3372291667    0.0000000000        1.#INF        1.#INF

λ echo 0 0 | bin\cct.exe -I +proj=latlong +pm=paris
  0.0000000000    0.0000000000        1.#INF        1.#INF # this should be  2.3372291667    0.0000000000

The forward operation is handled correctly because this

    case PJ_IO_UNITS_RADIANS:
        if (INPUT_UNITS==PJ_IO_UNITS_RADIANS)
            break;

happens in pj_fwd_finalize (line 131 and below), which basically omits the adjusted the longitude back to the greenwhich meridian. A similar construct is missing in pj_inv_prepare (remember we are going backwards now), so the longigude is adjusted to the desired prime meridan and then back to greenwich. The last bit breaks the test.

-------------------------------------------------------------------------------
tolerance 15 cm # lax tolerance due to widespread bad egm96 file
accept 12.5 55.5 0
expect 12.5 55.5 -36.0213

This comment has been minimized.

@kbevers

kbevers Jan 26, 2018

Member

With cs2cs I get the following output:

λ echo 55.5 12.5  0|bin\cs2cs.exe -f %.2f +proj=latlong +to +proj=latlong +geoidgrids=egm96_15.gtx  +axis=neu +ellps=GRS80
12.50   55.50 42.37

C:\dev\build\ninjaproj
λ echo 55.5 12.5  0|bin\cs2cs.exe -f %.2f +proj=latlong +geoidgrids=egm96_15.gtx  +axis=neu +ellps=GRS80 +to +proj=latlong
12.50   55.50 36.02

gie returns:

     -----
     FAILURE in 4D-API_cs2cs-style.gie(94):
     expected: 12.5 55.5 -36.0213
     got:      12.500000000   55.500000000   42.372592926
     deviation:  78393.893 mm,  expected:  150.000 mm

which corresponds to the first call to cs2cs above. I guess this all comes down to what we mean by +axis. Does +axis describe the ordering of the input or the output coordinates?

@kbevers

kbevers Jan 26, 2018

Member

With cs2cs I get the following output:

λ echo 55.5 12.5  0|bin\cs2cs.exe -f %.2f +proj=latlong +to +proj=latlong +geoidgrids=egm96_15.gtx  +axis=neu +ellps=GRS80
12.50   55.50 42.37

C:\dev\build\ninjaproj
λ echo 55.5 12.5  0|bin\cs2cs.exe -f %.2f +proj=latlong +geoidgrids=egm96_15.gtx  +axis=neu +ellps=GRS80 +to +proj=latlong
12.50   55.50 36.02

gie returns:

     -----
     FAILURE in 4D-API_cs2cs-style.gie(94):
     expected: 12.5 55.5 -36.0213
     got:      12.500000000   55.500000000   42.372592926
     deviation:  78393.893 mm,  expected:  150.000 mm

which corresponds to the first call to cs2cs above. I guess this all comes down to what we mean by +axis. Does +axis describe the ordering of the input or the output coordinates?

@kbevers

This comment has been minimized.

Show comment
Hide comment
@kbevers

kbevers Jan 29, 2018

Member

Seems like this is getting close to the finish line. I still get errors when running the 5201.gie test. I guess it is not run on CI yet. The problem is related to +proj=geocent, which at the moment is not taken care of in the cs2cs emulation. I don't think it is too difficult to take care of.

Member

kbevers commented Jan 29, 2018

Seems like this is getting close to the finish line. I still get errors when running the 5201.gie test. I guess it is not run on CI yet. The problem is related to +proj=geocent, which at the moment is not taken care of in the cs2cs emulation. I don't think it is too difficult to take care of.

busstoptaktik and others added some commits Jan 5, 2018

Introduce compatibility for cs2cs-style proj-strings into the 4D API.
Parameters such as towgs84, nadgrids and geoidgrids was previously only
handled by pj_transform(). This commit add a compatibility layer in
proj_create() by calling the pj_cs2cs_emulation_setup() function. This
function sets up a handful of predefined transformation objects on the
PJ object that is being created. Each of these transformation objects
are related to the cs2cs-style parameters we are trying to emulate in
the 4D API. That is, if the +towgs84 parameters is used we create
P->helmert with the parameters specified in +towgs84. Similarly for
+axis, +nadgrids and +geoidgrids. When these transformation objects
exists we use them in the prepare and finalize functions in pj_fwd/
pj_inv. If no cs2cs-style parametes are specified we skip those
parts of the prepare and finalizing steps.

Co-authored-by:Thomas Knudsen <thokn@sdfe.dk>
Co-authored-by:Kristian Evers <kristianevers@gmail.com>
Test material for the cs2cs emulation in the 4D API.
The GIGS tests that are known to work are added to the CMake test
setup. The GIGS gie files have been auto-translated from the
existing json-files and some corrections to tolerances have been
necessary since gie uses different norms than GIGS specify. The GIGS
tolerances are specified as the infinity norm of angular coordinates,
whereas gie uses the actual distances between calculated and expected
coordinates (using geodesics). In a few tests +towgs84 is overriden
from the EPSG inits to avoid creeping numerical inaccuracy in
roundtrips.

Co-authored-by: Thomas Knudsen <thokn@sdfe.dk>
Co-authored-by: Kristian Evers <kristianevers@gmail.com>

@busstoptaktik busstoptaktik merged commit 8f3345f into OSGeo:master Feb 1, 2018

3 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.3%) to 74.85%
Details

@busstoptaktik busstoptaktik deleted the busstoptaktik:4D-API_cs2cs-style branch Feb 1, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment