Plans to support 4D coordinates?

[kbevers]

Now that pyproj uses PROJ 6 behind the curtains it would be nice to be able to do more than just plain projections and simpel datum shifts. So, are there any plans to support 4D coordinates? Without having studied the pyproj in details it seems to be fairly straight forward to me to extend pyproj.transform() with an optional t input.

This would allow users of pyproj to do transformations that uses the 14-parameter version of the Helmert transform. Personally I would like to use it so that I can support the NKG transformations code based on pyproj.

[snowman2]

Sounds like a good idea to me. :+1:

[snowman2]

Mind adding an example here to be used for development purposes?

[kbevers]

Sounds like a good idea to me. 👍

Great!

Mind adding an example here to be used for development purposes?

Sure. Let's start with something basic that should be supported by PROJ out-of-the-box:

echo 3513638.19380   778956.45250    5248216.46900 2008.75| cct +init=ITRF2008:ITRF2000  3513638.1999    778956.4533  5248216.4535     2008.7500

[snowman2]

I am unfamiliar with 4d transforms, are there cases where there is a CRS to CRS transformation with time as well?

[snowman2]

Also, have an initial implementation in PR #218.

[kbevers]

I am unfamiliar with 4d transforms, are there cases where there is a CRS to CRS transformation with time as well?

Yes. There's a few examples in https://github.com/OSGeo/proj.4/issues/1354 which I've just reported while trying to come up with some test material for you. I think that issue is only related to cs2cs so hopefully you can use the two transformations from that issue for your tests.

Also, have an initial implementation in PR #218.

Awesome, thanks!

[snowman2]

I added it in here, however the transform using CRS to CRS does not seem to do anything (refer to the tests).

[snowman2]

From @kbevers:

This doesn't seem right. From PROJ master:

$ echo 3496737.2679 743254.4507 5264462.9620 2019.0 |./cs2cs EPSG:7789 EPSG:8401
3496737.76  743253.99 5264462.70 2019.0

[snowman2]

>>> transtime = Transformer.from_proj("EPSG:7789", "EPSG:8401")
>>> transtime.transform(3496737.2679, 743254.4507, 5264462.9620, 2019.0)
(3496737.2679, 743254.4507, 5264462.962, 2019.0)

[snowman2]

>>> pj = Proj("EPSG:7789")
>>> pj.crs
<CRS: epsg:7789>
Name: ITRF2014
Ellipsoid:
- semi_major_metre: 6378137.00
- semi_minor_metre: 6356752.31
- inverse_flattening: 298.26
Area of Use:
- name: World (by country)
- bounds: (-180.0, -90.0, 180.0, 90.0)
Prime Meridian:
- longitude: 0.0000
- unit_name: degree
- unit_conversion_factor: 0.01745329
Axis Info:
- Geocentric X[X] (geocentricX) EPSG:9001 (metre)
- Geocentric Y[Y] (geocentricY) EPSG:9001 (metre)
- Geocentric Z[Z] (geocentricZ) EPSG:9001 (metre)

>>> pj.crs.srs
'epsg:7789'

[snowman2]

The part that has me confused is this seems to work: https://github.com/pyproj4/pyproj/blob/623cbc8205c420e7f3a0f72a8f0410af430142a2/unittest/test_transformer.py#L114-L121

[snowman2]

@kbevers, one main difference I see is that pyproj is using proj_trans_generic and the code for cs2cs is using proj_trans.

[kbevers]

one main difference I see is that pyproj is using proj_trans_generic and the code for cs2cs is using proj_trans.

proj_trans_generic uses proj_trans internally so that is not enough to explain this :-) Maybe you should check the return code from proj_trans_generic to determine if you get back what you expect to get back.

[snowman2]

proj_trans_generic uses proj_trans internally so that is not enough to explain this :-)

Actually, it might. I did a test and added the proj_trans functionality to the Transformer class here, and everything works as expected.

>>> from pyproj import Transformer
>>> transformer = Transformer.from_proj("EPSG:7789", "EPSG:8401")
>>> transformer(3496737.2679,743254.4507,5264462.9620,2019.0)
(3496737.757717311, 743253.9940103051, 5264462.701132784, 2019.0

[snowman2]

Here is comparing both outputs:

>>> from pyproj import Transformer
>>> transformer = Transformer.from_proj("EPSG:7789", "EPSG:8401")
>>> transformer(3496737.2679,743254.4507,5264462.9620,2019.0)
(3496737.757717311, 743253.9940103051, 5264462.701132784, 2019.0)
>>> transformer.transform(3496737.2679,743254.4507,5264462.9620,2019.0)
(3496737.2679, 743254.4507, 5264462.962, 2019.0)

Both transformations use the same PJ*.

[snowman2]

Also, no error was raised by PROJ:

>>> transformer.transform(3496737.2679,743254.4507,5264462.9620,2019.0, errcheck=True)
(3496737.2679, 743254.4507, 5264462.962, 2019.0)

[kbevers]

Right, okay, that is odd. What does proj_trans_generic return? It should return the number of successfully transformed coordinates. There's a few cases where it can return 0 without setting an errno, maybe that is what is going on?

Is the same code used in the case of 3D transformations? Or is this a special code path in case of 4D coordinates?

[snowman2]

There's a few cases where it can return 0 without setting an errno

I didn't know that. I added a check for that in #228.

What does proj_trans_generic return?

It returns nothing. No error raised.

Is the same code used in the case of 3D transformations?

Yes, it is the same code (if you are referring to proj_trans_generic usage).

[kbevers]

It returns nothing. No error raised.

What does nothing mean? 0?

[snowman2]

Yeah, zero.

[kbevers]

There's two ways proj_trans_generic can return 0 without setting an error:

If either nullptr==P or 0==nx+ny+nz+nt is satisfied.

So I would check to see if that is the case.

[snowman2]

nullptr==P is not the case -> already checked with the proj_trans case.

[snowman2]

I think I found the issue, it is skipping because of equivalence logic. I will debug later.

[snowman2]

Yep, that was it. Fixed in #231. Thanks for the debug assistance!