Transformation pipelines #388
Conversation
…rmation-pipelines
Still lacking a lot, but the proof-of-concept has been accomplished
Much improved logic. Stuff takes shape!
The COORDINATE union eliminates explicit type-punning and makes the code much more clear
Each step is now introduced with the +step keyword. +first and +then are gone
The pipeline can now signal to the proj executable, that it provides angular output, and hence ask proj to avoid treating output as linear units (and hence multiply by the semimajor axis length). Proj.c now listens to what the pipeline manager says. This needs more work to support the -I (invproj) option: We need to specify both expected input and provided output, and to provide an override facility for the user.
…rmation-pipelines
As part of the pipeline work, reimplement the geodetic-to-geodetic functionality in proper 3D PJ style. The existing PJ_geodetic is mostly a stub, deferring the actual work to behind-the-scenes code in pj_transform.c. For datum shift pipelines, we want a proper 3D, accessible version. This is implemented in PJ_cart.c
The proj.c code is quite inpenetrable, so it takes quite some time to get things working. Now output is correct, but I assume that it can be made working in a less inelegant way...
Misc. Cleanups for pipeline. Debut for the 3D transformation program TRAN (still very rudimentary, but aims at becoming the 3D equivalent of PROJ. For now mostly a demo application for transformation pipelines.
License and descriptive headers added to the TRAN and transformation pipeline stuff
Added a handful of debugging tools for pipelines - and used them to eliminate a number of showstopping bugs. Also, incidentally, cleared up the flow of pj_fwd3d.c, which was unnecessarily convoluted and, additionally, appeared to have suffered under the encounter with a too agressive automated reformatting tool
pj_cart now again uses the code from geocent.c (which, having many years of use, is probably more robust than the new implementations). Also, finished the "add" pipeline functionality (which still lacks regression testing)
As part of the Helmert implementation, also added a rotations element (OPK: omega, phi, kappa) to the COORDINATE union (which might really better be called PJ_TRIPLET, although this may be slightly confusing since it also includes the two pairs XY and LP). Also returned PJ_cart.c to use the HM transformations (but eliminated the polar singularities): The pj_fwd/pj_inv i/o-mangling made the geocent.c stuff return wrong values.
Avoid hving arguments checked for functions that return same type as they accept.
And multiple minor improvements
No idea how it popped up again. Perhaps ctrl-y?
|
This would be material worth posting on the mailing list ! A few remarks and questions :
|
Tests was running on Windows, so PJ_pieline.h was found even when including pj_pipeline.h
|
Even @rouault , thank you for some very useful remarks - I will try to answer and/or comment. First things first... regarding pivot datum or not: when going from cartographic level accuracy to geodetic accuracy, it is actually rather the primary case not to have a pivot datum - the national realizations of ETRS89 are based on differing realizations to differing epochs, and you run into rather hairy operations in order to transform between those, even though they are all official ETRS89 realizations. Also (at least my experience from Danish, Faroese and Greenlandic systems), transformations between older systems can be in the form of polynomial corrections taking the older plane coordinate system to ED50/utm32 - which can then afterwards be transformed further to ETRS89 given access to the relevant Helmert parameters for the area in question (ED50 was rather jerky and transformation parameters are only sufficiently accurate for fairly small, typically national scale or less, areas). The use of the designation “WGS84” for the pivot datum (through the +towgs... keyword) is also somewhat unfortunate - because which of the WGS84 realizations, e.g. WGS84(G730, G873, G1150, G1674), is the pivot? It indicates that for cartographic purposes we may take ITRFxx=WGS84(Gxxxx)=ETRFxx=IGSxx - and at some time in history and for some applications that is certainly true. But 30 years after the first realization of WGS84, your average continent has moved by around 1 m, indicating that the identity between the realizations mentioned is now only true if you are happy with meter level accuracy. And for most civil uses (but probably not including road navigation), that is certainly not good enough. Hence, for greater accuracy, you need to know what you do. And while cs2cs works wonders in guessing the right incantation, I believe that as the accuracy requirements grow, we will need to provide a more tunable interface (not in order to take over the cs2cs and proj corners of the cartographic ecosystem, but rather to expand the range of validity and usefulness for libproj). Regarding PJ_cart.c and the question about geocentric etc: pj_geocent is a no_op - it just returns what was given to it and the actual transformation is carried out inside the stomach of pj_transform.c, using the routines in geocent.c, which can also be used from PJ_cart.c. I decided to implement my own routines as well, because I could find no solid references to the algorithm used in geocent.c, so I used the Bowring stuff, mentioned in the Hofmann-Wellenhof & Moritz reference in the code. It is well known and fairly simple to implement, so to have a “second opinion” from an algorithm I trusted, I implemented that one. But remember that this is work in progress, and I see no deep reason to keep my code in the codebase if I can find evidence that the geocent.c stuff is more trustworthy. A quick scan of Wikipedia does, however, indicate that both Bowring’s algorithm, and the one in geocent.c has been superseded by newer, faster, and more accurate algorithms - so perhaps it is worth reconsidering the contents of geocent.c As regards your comment that “geocentric” is a better name than “cartesian”, I tend to disagree: Cartesian indicates directly that the coordinates refer to a rectangular system, while geocentric just says something about the origin of the system: It could equally well be a system using geocentric latitude, longitude and height. Full semantic coverage is first reached by calling it “geocentric, cartesian”, which I’m sure we’ll both agree is way too long. I chose “cart” because in geodetic software, the indicators cart and crt are common - but to me “crt” still means cathode ray tube, and cart came with the added bonus that you can actually pronounce it as a word - so cart it was :-) About +proj=pipeline: It seems awkward to you at first, and it certainly also did to me: My first intention was to introduce a new keyword “+tran” for transformation (all projections are transformations, but not all transformations are projections), but since it would really just be a synonym for +proj anyway, I decided to defer it, thinking that the ugliness would provoke some good thoughts for improved syntax. That being said, the current implementation is technically very attractive, being simple and non-intrusive: Everything is done inside the PJ_pipeline.c file. I would prefer having at least a comparable degree of simplicity for any newly proposed syntax. About the public/private limit: I think we can greatly improve the public facing interface struct PJconsts; That way, for user code, PJ will be entirely opaque, not posing any risk of accidentally breaking anything - user code will only have to ferry around pointers to PJ-objects, so it does not need to know the internals, which will be declared in projects.h. Obviously, for compatibility, we will have to keep the projPJ typedef, although I personally find the name rather ugly, and believe we should try to stay in the PJ/pj namespace, except for the XYZ, LPZ, etc. data types, which are perfectly generic geodetic datatypes, which should not be considered private to libproj, but as a general service to the entire community, delivered through a lightweight proj_api.h But while this is related, it is also somewhat orthogonal to the pipeline stuff, where I agree that parts of PJ_pipeline.h should be moved to proj_api.h - and the rest to projects.h: I really only introduced PJ_pipeline.h to keep the experimental stuff somewhat localized. Your last comment about C-API versus “+proj=pipeline”: The code is prepared for cs2cs-style use. But my assumption is that the ordinary user does not know about the details of pipelines, but rather just uses an init-file provided by his national geodetic authority, such that e.g. the Danish transformation from the old System34 to ETRS89/utm32 can be invoked as “tran +init=DK:s34j_utm32etrs”, or some such. In my mental model, the C-API to pipelines is simply pj_init, and the pipeline syntax, and the pipeline primitives, make that API configurable. Does that sound reasonable to you? Thank you again for some very interesting and constructive comments /thomas |
Added pj_is_isomorphic to proj.def
Repaired the inverse transformation, added "transpose" and "approx" options, and described their effect in a long comment elaborating on the mechanics of the rotation matrix, and the sign conventions in common use.
... and remove accidentally added file PJ_helmert.c, in preparation of renaming PJ_pipeline_eæements.c to PJ_helmert.c
Added new src/makefile.vc from Kristian Evers @kbevers, updated build system files to reflect PJ_pipeline_elements.c -> PJ_helmert.c rename
... and a few improved comments in projects.h
Apart form being the generic header for the TRIPLET union, which is useful for geodetic applications in general, proj.h also provides access to the minimalistic API consisting of: pj_init_plus pj_apply pj_free i.e. the bare essentials needed for using proj.4 in a TRIPLET based workflow. I plan to add a few things to this, while keeping it both simple and useful
(adding amissing .c)
In MSVC HUGE_VAL is not a compile time constant, so we cannot have a constant initialized to HUGE_VAL. Intead introduce a generic TRIPLET constructor, pj_triplet.
(not ready for prime time)
…rmation-pipelines
|
Norm Olsen describes a related effort for the CS-Map library over at https://trac.osgeo.org/csmap/wiki/CsMapRfc2 - worth reading! |
|
Hi @busstoptaktik, that is very interesting work you're doing! I've been looking into ways to define a coordinate system which is Gauß-Krüger Zone 4 plus an affine 3D transform (to define an arbitrarily oriented local tangent space) and found this branch. I compiled proj, recompiled stuff that depends on it (e.g., PostGIS) and played with the "tran" tool and the Helmert transform. While I could finally get +proj=helmert to apply my affine transform (which is defined by a 2x2 rotation matrix and a translation vector, so it was necessary to compute the angle rz, convert the angle to arcseconds, as well as deriving the proper scale), I would really love to see a generic 3D affine transform that allows specifying the six unique elements of a 3x3 rotation matrix R00, R01, R02, R11, R12, R22 and the 3-vector for translation. That would allow you to easily concatenate multiple transformations in local Cartesian spaces until a final transform like tmerc describes how the datum is transformed to WGS-84. So having "tran" work, I was curious if this is going to work in PostGIS. I defined a custom SRID with a pipeline to see if I can add an extra Helmert transformation on top of GK-Z4. So I define a spatial reference system and try to use it with ST_Transform: INSERT INTO spatial_ref_sys VALUES(1000, 'test', 1000, NULL, '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs');
SELECT ST_AsText( ST_Transform( ST_GeomFromText('POINT(4470814.627454 5640426.090921)', 1000), 4326) );
>>> "POINT(11.5836498596428 50.898421786718)"(not exactly what I expected: 11.585083007813, 50.899658203125, but close enough for now...) As soon as I turn that into a pipeline (even with only a single transform), the result is very different: INSERT INTO spatial_ref_sys VALUES(1000, 'test', 1000, NULL, '+proj=pipeline +step +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs');
SELECT ST_AsText( ST_Transform( ST_GeomFromText('POINT(4470814.627454 5640426.090921)', 1000), 4326) );
>>> "POINT(-13.4108026989403 0.00376315356350652)"Do you have any hints for me if what I want to achieve is at all possible with your modifications and how to get this to work? |
|
@rhuitl Thanks for your comments. And thanks for testing the transformation pipelines in the real world. We haven't even reached that point yet. Please continue to test the code when new commits appear. Bug reports and suggestions on how to improve the interface etc. are much appreciated! Your idea about using the rotation matrix directly in the pipeline is definitely worth considering. |
|
@rhuitl: Thanks for stress testing this attempt at proj.4 enhancement: I really appreciate that, as this PR is essential for implementing transformations for historical Danish, Faroese, and Greenlandic systems. Unfortunately I am busy with other tasks for the moment, and will not be back on this project for the next few weeks. I will do my best to follow up some time before mid-August. /Thomas |
|
@busstoptaktik That's no problem, I'll keep watching this PR and can redo my tests when you think it makes sense :) |
|
Do you guys have some cycles to polish this up and push it in to 4.9.3, or do you want to wait? |
|
It has to wait: It is actually quite intrusive, and needs extensive testing I do not want to delay 4.9.3 by attempting to get this merged. I will, however, continue the work in a few weeks time. /Thomas Den 17/08/2016 16.11 skrev "Howard Butler" notifications@github.com:
|
|
Closing this PR - refactoring and resubmitting as a series of smaller PRs |
|
What is the status on implementing S34s/j in Proj4? Since TrLib is deprecated it would be nice to move to Proj.4. |
|
@haugwitz - yes, we have a working implementation of s34 (actually of s45, the s34 for Bornholm, but the code is identical to s34j/s, only the constants differ. I'm fully booked for the next two weeks, but please contact me directly at thokn AT sdfe DOT dk for more info about testing. |
With the completion of the moving of projection specific parameters to opaque objects, the size of the PJ object is now constant and deterministic.
This means that pj_init, which returns a pointer to the PJ object may, syntactically speaking, just as well return a pointer to an array of PJ objects.
This means that we may implement transformations as pipelines of steps, just by introducing one new projection, which handles the handover between a series of operations, each implemented as individual projections.
That is the approach taken in this pull request, which introduces a new 3D transformation program “tran”, the “pipeline” projection handler, the “step” keyword, and reinroduces the “inv” qualifier.
It also introduces a number of transformation primitives (each implemented as a PJ-style projection):
+proj=cart: For geographic / 3D cartesian transformation
+proj=helmert: For 7 parameter Helmert shifts
+proj=add: For 3 parameter shifts
(more will follow)
And implements a 3D transformation program "tran". Details, User’s Guide and binary downloads are available over at http://thomasknudsen.net/proj4.html
The pull request also introduces a more semantically meaningful approach to type-punning than the “use preprocessor redefinitions to make both LPZ and XYZ look like an UVW” approach used currently: The new approach introduces the COORDINATE union (which should perhaps be renamed TRIPLET or PJ_TRIPLET), which allows us to have all the flexibility of the UVW-approach, while still keeping the semantic meaning which is so much more evident when the data type actually reflects the data contained.
This is work in progress, only 10 days into its existence. The pull request is certainly not ready for merging at this stage, but there is a number of things I would like to discuss, and find a PR the most straightforward way of starting that discussion.
Discussion points:
Is this at all interesting? In my book, it introduces some very useful functionality in a almost entirely non-intrusive way. But different people, different bananas...
ABI breakage: It has become clear that a number of additions to the PJ data type would be beneficial. Would it be OK to introduce these in 4.10.0? Currently I have implemented workarounds, but added elements would simplify the code, and open up new possibilities.
If this is interesting, what about the syntax? Any advice for improvements? One thing could be to make +step= a synonym for +proj, so the syntax could be simplified from “+step +proj=foo...” to “+step=foo...”
4.
Your concern here...