Do we really support gridshifts on cartesian input?

Not in "classic PROJ" mode. This can only occur if when using the deformation operation which is not compatible with classic PROJ.

So this will first really be an issue to decide on when we introduce global +inv handling and/or a proj=cartesian (or proj=geocent in deprecated PROJ terms).

I'll leave the code in, but we will have to decide what we mean with this. I can imagine situations where it makes sense, but mostly for convenience (but then again, there is nothing in this PR which cannot be handled by properly constructed pipelines - it is by and large a matter of improving backward compatibility by providing some convenience functionality, for interoperability between new code and old parameter file systems)

(but then again, there is nothing in this PR which cannot be handled by properly constructed pipelines - it is by and large a matter of improving backward compatibility by providing some convenience functionality, for interoperability between new code and old parameter file systems)

That's my point. Keep it as simple as possible and don't open the door too much for new ideas based on the classic PROJ.4 idioms.

OK, I get your point now, and I fully agree, so I will remove that chunk of code

pj_param_exists?

Yes, it makes sense here, and we could redesign the entire parameter handling in this function to use pj_param_exists, which I think we should do, but not now. Changed this instance, because it is different in having its work already done by a friend up in pj_datum_set

in pj_inv you define INPUT_UNITS and OUTPUT_UNITS. Doing that here as well will make it easier to read the code

Yes, and I will do that. I managed to do the mental math necessary in the forward case, without having this translation/reminder, whereas when having to do the same thing backwards and in high heels, in the inverse case, I really needed these persistent pings

Use PJ_INV here. And many other places in the code below.

I don't agree. PJ_INV/PJ_FWD were introduced for your convenience. I actually find it quite a lot more readable to have negative mean backward, and positive mean forward.

And since there is more than enough of layered negations in this code, I really prefer the concise to the verbose here.

I guess we will just have to agree to disagree here :-)

Sounds very agreeable :-)

That being said, if you find it unreadable, I would not mind switching to PJ_FWD/PJ_INV once we have established that the material works as intended (last step before merging).

It's just that I find the PJ_FWD/PJ_INV convention less readable than the mathematical convention of +1/-1, and this stuff is pretty obnoxious already.

PJ_FWD here and below.

(ibid.)

At this point I believe P->axis can be expected to contain valid data if it is set. So no need for the strcmp.

The strcmp returns true (i.e. non-null) if *(p->param) != "enu", so it's a way to avoid calling axisswap if the data are already in the expected order.

I agree that the data must be valid at this point, so this is not a validation, just a cheap optimization, which may deserve a comment for explanation.

Ah of course. Makes sense.

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.

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.

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.

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.

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 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).

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

Duly noted.

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.

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

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

May I suggest calling it coord instead. I think it reads better.

coo has been conventionally used here, there, and everywhere:

grep coo *.c *.h | grep -v coord | wc -l

gives 74 hits.

coord is used in methods and type names (PJ_COORD, proj_coord (), etc.) - coo is used for an otherwise unspecific instance, intentionally modelled after the equivalent unspecificness of the versatile variable foo of foo, bar, baz fame.

coord is used as an instance name a few places in proj_4D_api.c and gie.c - I would be less reluctant renaming these instances to coo, than changing the already widespread coo idiom.

coo has been conventionally used here, there, and everywhere

I know, but I don't like the convention. It somehow made sense previously when coo was also used in the now gone PJ_OBS where a short name was justified. Since that no longer exists I think it is perfectly reasonable to use a more descriptive convention.

I would also object to using non-descriptive names like foo and bar in anything other than example code so I am not buying that argument.

Anyway, my point is just that coord makes more sense to readers unfamiliar with the code-base as a whole. I am not suggesting that the existing occurences of coo be replaced now, I just wanted to raise a flag here since it has been bothering me for a while. Do what you feel is best - ultimately it is in the nitpicking department!

@kbevers: Since you dislike the convention, I think we should consider reconsidering the convention(s).

From classic PROJ/Gerald, we have P for the generic PJ object. From 3D-PROJ/Frank, we have ctx for the generic PJ_CONTEXT object. So apparently, the smaller the object, the longer the name.

With that in mind, it actually makes quite a lot of sense to switch coo (same length as ctx, but smaller) to your preferred coord (longer and smaller than ctx).

But if that's to be done I suggest we do it consequently, and reserve a rainy summer day, in time for 5.1.0, for a sweep-over of the entire code base, switching coo to coord - and collecting a list of other paper cuts/minor annoyances to be remedied at the same occasion.

I'll buy the beer, and open the issue for discussion - you bring the stamina, and I'll try to follow your typing pace :-)

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?