search

OSGeo / PROJ

PROJ - Cartographic Projections and Coordinate Transformations Library
https://proj.org
Other
1.75k stars 791 forks source link

New API applies inverse (wrong) +towgs84 helmert transform #783

Closed c0nk closed 6 years ago

c0nk commented 6 years ago

Consider the transform from WGS84 to CH1903+/LV95. Building the transform with proj_create_crs_to_crs(ctx, "epsg:4326", "epsg:2056", 0); yields wrong results. Apparently helmert transform is being applied in the wrong direction. Alternatively, if I create the transform using the epsg definition it has the same issue:

proj_create(ctx,
    "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 "
    "+x_0=2600000 +y_0=1200000 +ellps=bessel "
    "+towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs");

However, I can get the correct transform with the new API by negating the +towgs84 params:

proj_create(ctx,
    "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 "
    "+x_0=2600000 +y_0=1200000 +ellps=bessel "
    "+towgs84=-674.374,-15.056,-405.346,0,0,0,0 +units=m +no_defs");
busstoptaktik commented 6 years ago

Hello @c0nk - I have actually been working on this tonight, and will be wrapping the material up in a number of pull requests during the next few days.

The intention is to make the proj_create_crs_to_crs function work as expected, i.e. to return a PJ object which, when used in connection with the proj_trans family of functions, make them work as a corresponding call to the old pj_transform would.

The way you got things working, using proj_create is possible, but not recommendable: The thing is that while the old API is based on a declarative programming style (i.e. describing two systems, and leave it to PROJ to figure out how to transform between them), the new API is more in the style of imperative programming (i.e. describe the steps needed to go from system A to system B, without necessarily mentioning any of these systems explicitly).

In the long run, the intention with the proj_create_crs_to_crs function is to bridge the gap between the declarative and the imperative by, behind the scenes, looking up imperative recipes for going from system A to system B, without needing the WGS84 pivot datum, which in PROJ (geodetically speaking) is ill defined.

So when using proj_create you should write up a pipeline, daisy chaining the steps you need to take to go from the input system to the output system.

Pipelines tend to get long and messy (they execute pretty fast, though), and their primary use is to let geodetic authorities, such as National Mapping Agencies, distribute transformations in a reasonably human readable way, by use of PROJ parameter files.

But thanks for your active participation in this release cycle - I hope you are willing to follow up when I have wrapped up and merged tonight's work (probably some time during the weekend)

c0nk commented 6 years ago

@busstoptaktik Thanks for the quick reply and the good work. Why is my use of proj_create not recommended? It is my understanding that the EPSG definitions use WGS84 as a "pivot". So it should be reasonable to use the EPSG definition directly if the source or destination CRS is WGS84. Or am I missing something?

By including +towgs84 in the definition string, it is basically saying "use WGS84 as the source CRS and apply the helmert transform as a preprocessing step before applying the primary transformation". Looking at it this way, I expect

proj_create(ctx,
    "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 "
    "+x_0=2600000 +y_0=1200000 +ellps=bessel "
    "+towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs");

to also give the correct result. Can you verify this will be the case after your patches get merged?

c0nk commented 6 years ago

I expect the following two transformations to be equivalent:

proj_create(ctx,
    "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 "
    "+x_0=2600000 +y_0=1200000 +ellps=bessel "
    "+towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs");

and

proj_create(ctx,
    "+proj=pipeline "
    "+step +proj=cart +ellps=WGS84 +no_defs "
    "+step +proj=helmert x=674.37400 y=15.05600 z=405.34600 +no_defs +inv "
    "+step +proj=cart +ellps=bessel +no_defs +inv "
    "+step +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 "
    "+k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +units=m +no_defs ");
busstoptaktik commented 6 years ago

Can you verify this will be the case after your patches get merged?

I cannot verify anything right now (need some sleep). Verification when I wrap up the work will be easier if you, in the meantime, post a few points with input coordinates and expected output coordinates.

c0nk commented 6 years ago

7.438632495 E, 46.951082877 N -> x=2600000.0, y=1200000.0 (sub meter accuracy).

c0nk commented 6 years ago

Just let me know when you have something and I will integrate / verify your new changes as well.

busstoptaktik commented 6 years ago

@c0nk I have now pushed the first version of the work - it still does not include the "working around proj_def.dat" thing, which may come in a separate PR later today or tomorrow.

Before pushing on, I need to see the results from the CI builds - and would also be happy to see any comments or observations you may have at this stage.

busstoptaktik commented 6 years ago

@c0nk, I have now checked your test material by putting it in a file c0nk.gie:

<gie>
-------------------------------------------------------------------------------
operation  proj=somerc
           lat_0=46.95240555555556 lon_0=7.439583333333333 k_0=1
           x_0=2600000 y_0=1200000 ellps=bessel
           towgs84=674.374,15.056,405.346
-------------------------------------------------------------------------------
tolerance 20 cm
accept    7.438632495 46.951082877
expect    2600000.0  1200000.0
-------------------------------------------------------------------------------
</gie>

Then running:

gie c0nk.gie

With this result:

-------------------------------------------------------------------------------
Reading file '..\proj.4\test\gie\c0nk.gie'
-------------------------------------------------------------------------------
total:  1 tests succeeded,  0 tests skipped,  0 tests failed.
-------------------------------------------------------------------------------

If you do not mind, I would like to incorporate this test into the 4D-API_cs2cs-style.gie set.

c0nk commented 6 years ago

@busstoptaktik I haven't got a chance to verify your work yet because I am getting a compile error with the latest code. See #781 .

busstoptaktik commented 6 years ago

@c0nk I have now merged #791, so you can use proj_context_errno and do not need to include proj_api.h

c0nk commented 6 years ago

The following is working now as expected:

proj_create(ctx,
    "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 "
    "+x_0=2600000 +y_0=1200000 +ellps=bessel "
    "+towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs");

However, the pipeline variant:

proj_create(ctx,
    "+proj=pipeline "
    "+step +proj=cart +ellps=WGS84 +no_defs "
    "+step +proj=helmert x=674.37400 y=15.05600 z=405.34600 +no_defs +inv "
    "+step +proj=cart +ellps=bessel +no_defs +inv "
    "+step +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 "
    "+k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +units=m +no_defs ");

is now failing with error major axis or radius = 0 or not given.

kbevers commented 6 years ago

The problem seems to be +proj=helmert that suddenly requires a ellipsoid to be set. If you remove +no_defs from the helmert step you can proceed. Although there's still a bug to fix.

busstoptaktik commented 6 years ago

Probably a logical bug in my update of pj_init yesterday. Will look into it later today. Basically we should probably do as in the pipeline case and force WGS84 if no_defs set

busstoptaktik commented 6 years ago

Probably a logical bug in my update of pj_init yesterday. Will look into it later today. Basically we should probably do as in the pipeline case and force WGS84 if no_defs set

kbevers commented 6 years ago

It's a problem for all operations that has needs_ellps == 0 i think.

c0nk commented 6 years ago

If an operation doesn't need an ellipse. Why does it need to default to WGS84?

kbevers commented 6 years ago

If an operation doesn't need an ellipse. Why does it need to default to WGS84?

Because we've bolted a bunch of new stuff on top of an old architecture where it made sense to have such a requirement. We are gradually refactoring but it is difficult as long as we still have to respect old APIs. So once in a while something is overlooked like in this case.

busstoptaktik commented 6 years ago

I think I (finally) nailed it - also for pipelines - in PR #794

c0nk commented 6 years ago

I have tested with a557b49e85ebfbfe0af076bc529002ca60e4a045. My pipeline example from above is still not working for me. Same error.

busstoptaktik commented 6 years ago

The correction is in #794 (and partially in #793)

c0nk commented 6 years ago

Ah OK. I will test that when it merges.

busstoptaktik commented 6 years ago

It may take some time: The CI queues are apparently heavily loaded

busstoptaktik commented 6 years ago

Apparently not entirely done - still some problems when PROJ_LIB not set.

c0nk commented 6 years ago

I have tested on master and the original issues are solved now.

c0nk commented 6 years ago

I think I have reproduced your PROJ_LIB issue with +towgs84.

busstoptaktik commented 6 years ago

Yes - I believe it is solved now, in 715ec0b90aec05f15554c6591ea8cfd09d5cf043 (checked using your test example). Can you confirm that?