Discussion:
[Proj] Understanding 2D Helmert
Andre Joost
2018-10-26 13:43:03 UTC
Permalink
Dear list,

since PROJ 5.2 now has made its way into the OSGEO4W binaries, I try to
understand the new features.

I started up with the 2D helmert transformation with this command line:

echo 1000 0 | proj +proj=helmert +convention=coordinate_frame +x=0 +y=0
+s=1. +theta=-5

This should rotate a point with X=1000 Y=0 meters by theta=-5 around the
origin.

But the result of "17.45 0" does not really make sense to me.
X=2000 results to a new X=34.91.

Can anybody shed a light on this?
Tobias Wendorff
2018-10-26 13:52:22 UTC
Permalink
Post by Andre Joost
But the result of "17.45 0" does not really make sense to me.
X=2000 results to a new X=34.91.
Can anybody shed a light on this?
Seems to be in radians: 17.45 to degrees* = 999.8113525
*to degrees: 360/(2*PI()) * 17.45
Kristian Evers
2018-10-26 14:40:17 UTC
Permalink
If that is the case it should just be a matter of changing the
type of the output units internally in the operation. If you
submit a new ticket on the github page with a short description
I'll take a look at it within the next couple of days.

Thanks for reporting this.

/Kristian

-----Oprindelig meddelelse-----
Fra: proj-***@lists.maptools.org <proj-***@lists.maptools.org> På vegne af Tobias Wendorff
Sendt: 26. oktober 2018 15:52
Til: PROJ.4 and general Projections Discussions <***@lists.maptools.org>
Emne: Re: [Proj] Understanding 2D Helmert
Post by Andre Joost
But the result of "17.45 0" does not really make sense to me.
X=2000 results to a new X=34.91.
Can anybody shed a light on this?
Seems to be in radians: 17.45 to degrees* = 999.8113525
*to degrees: 360/(2*PI()) * 17.45
Even Rouault
2018-10-26 15:37:32 UTC
Permalink
Post by Kristian Evers
If that is the case it should just be a matter of changing the
type of the output units internally in the operation. If you
submit a new ticket on the github page with a short description
I'll take a look at it within the next couple of days.
Things are fine. This is just that using 'proj' for such operation is
inapproprate. You should use the new 'cct' binary.

'proj' is just for plain-old projection methods that take long, lat as input,
whereas here input coordinates are in the projected space

echo 1000 0 0 | cct +proj=helmert +convention=coordinate_frame +theta=-162000

outputs
707.1068 707.1068 0.0000 inf

(162000 arc-second = 45 deg)

(due to a current oddity of cct, you need to specify 3D coordinates, or add "-
z 0" in the parameters )
--
Spatialys - Geospatial professional services
http://www.spatialys.com
Kristian Evers
2018-10-26 17:47:03 UTC
Permalink
Ah yes, I missed that proj was used in the example.
The problem is the same, just located in proj.c instead.

I have pushed a fix here: https://github.com/OSGeo/proj.4/pull/1162

While proj is mainly meant to be used to project geodetic
coordinates it doesn’t hurt to treat other input correctly
when it is possible. It is a bit closer to doing that now. I
am sure there are some exceptions still that I haven’t
found (the code in proj.c structured quite poorly)

/Kristian

On 26 Oct 2018, at 17:37, Even Rouault <***@spatialys.com<mailto:***@spatialys.com>> wrote:

On vendredi 26 octobre 2018 14:40:17 CEST Kristian Evers wrote:
If that is the case it should just be a matter of changing the
type of the output units internally in the operation. If you
submit a new ticket on the github page with a short description
I'll take a look at it within the next couple of days.

Things are fine. This is just that using 'proj' for such operation is
inapproprate. You should use the new 'cct' binary.

'proj' is just for plain-old projection methods that take long, lat as input,
whereas here input coordinates are in the projected space

echo 1000 0 0 | cct +proj=helmert +convention=coordinate_frame +theta=-162000

outputs
707.1068 707.1068 0.0000 inf

(162000 arc-second = 45 deg)

(due to a current oddity of cct, you need to specify 3D coordinates, or add "-
z 0" in the parameters )
--
Spatialys - Geospatial professional services
http://www.spatialys.com
Andre Joost
2018-10-26 17:54:22 UTC
Permalink
Post by Even Rouault
Post by Kristian Evers
If that is the case it should just be a matter of changing the
type of the output units internally in the operation. If you
submit a new ticket on the github page with a short description
I'll take a look at it within the next couple of days.
Things are fine. This is just that using 'proj' for such operation is
inapproprate. You should use the new 'cct' binary.
Thanks, that worked for me in this way:

cct -z 0 +proj=helmert +convention=coordinate_frame +x=239228.435
+y=6693404.951 +s=1.0003267 +theta=18795.4443 <HelmertIn.txt
Post by Even Rouault
Post by Kristian Evers
HelmertOut.txt
Maybe this example should be added to the Helmert and/or cct man pages.

Greetings,
Andre Joost
Even Rouault
2018-10-26 18:11:03 UTC
Permalink
Post by Andre Joost
Post by Even Rouault
Post by Kristian Evers
If that is the case it should just be a matter of changing the
type of the output units internally in the operation. If you
submit a new ticket on the github page with a short description
I'll take a look at it within the next couple of days.
Things are fine. This is just that using 'proj' for such operation is
inapproprate. You should use the new 'cct' binary.
cct -z 0 +proj=helmert +convention=coordinate_frame +x=239228.435
+y=6693404.951 +s=1.0003267 +theta=18795.4443 <HelmertIn.txt
Post by Even Rouault
Post by Kristian Evers
HelmertOut.txt
Maybe this example should be added to the Helmert and/or cct man pages.
Where do those figures come from ? Is it some official transformation
registered somewhere ?
--
Spatialys - Geospatial professional services
http://www.spatialys.com
Andre Joost
2018-10-27 08:00:17 UTC
Permalink
Post by Even Rouault
Post by Andre Joost
cct -z 0 +proj=helmert +convention=coordinate_frame +x=239228.435
+y=6693404.951 +s=1.0003267 +theta=18795.4443 <HelmertIn.txt
Post by Kristian Evers
HelmertOut.txt
Maybe this example should be added to the Helmert and/or cct man pages.
Where do those figures come from ? Is it some official transformation
registered somewhere ?
See
https://gis.stackexchange.com/questions/300219/helmert-transformation
and
http://www.maanmittauslaitos.fi/sites/default/files/Finnish_Coordinate_Systems.pdf,
chapter 6.2

Seems they have 2D helmert transformations defined between KKJ and
ETRS-TM35 within every triangulation triangle. The link given in the
document is not valid anymore, but now available on
http://coordtrans.fgi.fi/kkj_EUREF-FIN.jsp after login.

One other thing that strikes me:
I have 8 lines within my input file, but only the first 7 get
transformed. The last input line does not end with CR LF.

No difference if I use
cct -z 0 -o cctout.txt +proj=helmert +convention=coordinate_frame
+x=239228.435 +y=6693404.951 +s=1.0003267 +theta=18795.4443 Helmertin.txt

cs2cs behaves correctly in this respect, all 8 lines transformed.
So I have to end my file with a blank line to get all points with both
apps. Tested only on Windows 7.

Greetings,
Andre Joost
Even Rouault
2018-10-27 22:03:57 UTC
Permalink
Post by Andre Joost
Post by Even Rouault
Post by Andre Joost
cct -z 0 +proj=helmert +convention=coordinate_frame +x=239228.435
+y=6693404.951 +s=1.0003267 +theta=18795.4443 <HelmertIn.txt
Post by Kristian Evers
HelmertOut.txt
Maybe this example should be added to the Helmert and/or cct man pages.
Where do those figures come from ? Is it some official transformation
registered somewhere ?
See
https://gis.stackexchange.com/questions/300219/helmert-transformation
and
http://www.maanmittauslaitos.fi/sites/default/files/Finnish_Coordinate_Syste
ms.pdf, chapter 6.2
Seems they have 2D helmert transformations defined between KKJ and
ETRS-TM35 within every triangulation triangle. The link given in the
document is not valid anymore, but now available on
http://coordtrans.fgi.fi/kkj_EUREF-FIN.jsp after login.
Thanks for the context. I see the EPSG dataset has a transform from KKJ
Geographic to ETRS89 using the parameters given in "6.1 The 3D transformation"

As far as the 2D helmert per triangle, this would be impractical to
incorporate that in the to-be PROJ database. But I guess someone could create
a grid that approximates the result of such transformation. Depending on the
step chosen, you could probably achieve very similar results to the "exact"
method.

One thing that strikes me in the description of the method is that you might
have some discontinuities when considering points on each side of an edge of
the triangulation. But perhaps the method is a bit more subtle than this
explanation and you don't just take the delta_x, delta_y, a1, a2, b1, b2
parameter of the triangle, but rather consider values of those parameters at
the 3 points of the triangle, and do a barycentric interpolation with the
coordinates of the point to be transformed. That way you would have
continuity.
Post by Andre Joost
I have 8 lines within my input file, but only the first 7 get
transformed. The last input line does not end with CR LF.
Would be worth opening a ticket about that.

Even
--
Spatialys - Geospatial professional services
http://www.spatialys.com
Kristian Evers
2018-10-28 15:02:05 UTC
Permalink
On 28 Oct 2018, at 00:03, Even Rouault <***@spatialys.com<mailto:***@spatialys.com>> wrote:

On samedi 27 octobre 2018 10:00:17 CEST Andre Joost wrote:
Am 26.10.18 um 20:11 schrieb Even Rouault:
On vendredi 26 octobre 2018 19:54:22 CEST Andre Joost wrote:
Thanks, that worked for me in this way:

cct -z 0 +proj=helmert +convention=coordinate_frame +x=239228.435
+y=6693404.951 +s=1.0003267 +theta=18795.4443 <HelmertIn.txt

HelmertOut.txt

Maybe this example should be added to the Helmert and/or cct man pages.

Where do those figures come from ? Is it some official transformation
registered somewhere ?

See
https://gis.stackexchange.com/questions/300219/helmert-transformation
and
http://www.maanmittauslaitos.fi/sites/default/files/Finnish_Coordinate_Syste
ms.pdf, chapter 6.2

Seems they have 2D helmert transformations defined between KKJ and
ETRS-TM35 within every triangulation triangle. The link given in the
document is not valid anymore, but now available on
http://coordtrans.fgi.fi/kkj_EUREF-FIN.jsp after login.


Thanks for the context. I see the EPSG dataset has a transform from KKJ
Geographic to ETRS89 using the parameters given in "6.1 The 3D transformation"

As far as the 2D helmert per triangle, this would be impractical to
incorporate that in the to-be PROJ database. But I guess someone could create
a grid that approximates the result of such transformation. Depending on the
step chosen, you could probably achieve very similar results to the "exact"
method.


Why is that? To me it seems somewhat similar to polynomial mappings (Horner),
in the way that many coefficients are likely to be stored for a given transformation.
Given a clever way of defining the triangles and their corresponding Helmet
parameters I think it should be possible to make these transformations work. At
least as PROJ strings.

One thing that strikes me in the description of the method is that you might
have some discontinuities when considering points on each side of an edge of
the triangulation. But perhaps the method is a bit more subtle than this
explanation and you don't just take the delta_x, delta_y, a1, a2, b1, b2
parameter of the triangle, but rather consider values of those parameters at
the 3 points of the triangle, and do a barycentric interpolation with the
coordinates of the point to be transformed. That way you would have
continuity.

I believe that is the case, yes. Last year we had a workshop within the
Nordic Geodetic Commision on the prospects of implementing the various
transformations used in the Nordic countries. This transformation was
brought up by both the Norwegians and the Fins and as far as I remember
it works exactly like you just described.
I know that both the Fins and the Norwegians are interested in seeing this
operations added in PROJ. I am sure they will be able to provide us with both
knowledge and test material should we want to implement this.

/Kristian


One other thing that strikes me:
I have 8 lines within my input file, but only the first 7 get
transformed. The last input line does not end with CR LF.

Would be worth opening a ticket about that.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com<http://www.spatialys.com/>
Even Rouault
2018-10-28 16:46:22 UTC
Permalink
Post by Kristian Evers
Why is that? To me it seems somewhat similar to polynomial mappings
(Horner), in the way that many coefficients are likely to be stored for a
given transformation. Given a clever way of defining the triangles and
their corresponding Helmet parameters I think it should be possible to make
these transformations work. At least as PROJ strings.
You would need to store hundreds of triangles and their parameters, and have
PROJ figure out in which triangle a point it. Certainly doable, but not
directly in the scope of my current work. And avoid the database to grow too
big with very particular transforms. That said, 600 triangles with 6 double
parameters each + the x,y coordinates is just 38.4 KB if you store that as a
compact binary blob. Some balance to know if it is something acceptable in-db
(EPSG style transformations have at most 18 parameters each) or out-db
(grids). The suggestion for the grid approach was to make it more immediatly
usable.
--
Spatialys - Geospatial professional services
http://www.spatialys.com
Luke Bergmann
2018-10-28 18:54:44 UTC
Permalink
Dear all,

As a fan of the proj[.4] project, as a lurker on these lists, and as hopefully a future contributor to proj, I thought I should jump in here with my first post. Please forgive me my mistakes in protocol and my ignorance! This thread's discussion has come close to functionality that a colleague of mine (David O’Sullivan) and I have long thought would really transform the capabilities of proj and projects downstream (e.g., QGIS), making them much more useful to geography, to ‘analytical cartography’, and to the digital humanities, among other pursuits.

In particular, we would like to suggest (and we are willing to contribute work to see that) proj implement transformations that are data-driven, or as Waldo Tobler called them, ‘empirical projections’. The easiest example would be the transformation associated with a cartogram. Of course, the math/heuristic to calculate a cartogram is clearly beyond the scope of a library like proj. However, each existing cartogram embodies a transformation, and having the ability to represent that transformation in an approximate/empirical/data-driven form within a proj string would enable your average GIS to suddenly work with a vastly expanded world of coordinate systems. Similarly, there is an explosion of creative transformations out there these days (think about all the heuristics unfolding the globe these days) that may themselves be data-driven and/or may not have been implementable in proj for various reasons—this would enable everyday GIS users to use them.

We recognize there is already a ‘grid’ functionality for datum transformations within proj, but having a proj projection that interprets and implements arbitrary projections encoded within proj strings is both closer to the spirit of understanding these transformations as projections (as Tobler and others have advocated) and would seem to be important from a practical perspective, as there are already ways in GIS such as QGIS to input proj strings and then call proj. As academics, we want this functionality to be usable by students in the introductory GIS course. There are other benefits to knowledge as well. The humanities, which see space as more complex, often imagined, would benefit from these transformations. And they would be key to proj enabling the ‘analytical' use of cartography that Waldo Tobler advocated with his ‘analytical cartography’.

How might this work? We gave a paper at SIGSPATIAL in 2017 that started to sketch out one idea. I’ve linked the paper below—I can email you a copy but the proj list doesn’t accept big attachments and I haven’t been able to set up the open sharing function with ACM yet, sorry! Note that the paper goes beyond what I’m able to describe here. Just as this list has been discussing transformations with hundreds of triangles, so too were we. To start, we thought that a proj string could encode how a set of points in one space map into coordinates within another space. Here would be a generic notional proj.4 CRS string:

+proj=ptobler +interpolation=interptype +mapping=mapstring

where:

- ptobler is a string literal that represents a possible shorthand name for the class of generic projections (‘Ptobler’ is a somewhat droll honorary offered to Waldo Tobler’s cartographic ideas which also references the cartographic achievements of Ptolemy. It first emerged out of conversations around the early Michigan Inter-University Community of Mathematical Geographers (MICMOG))

- interptype is a string representing the name of an interpolation method, e.g., polynomial, bicubic, thin plate spline, TIN. The interpolation method is used to estimate how arbitrary points in (φ,λ) coordinates map into (x,y) space.

- mapstring is a string serialization of the list of known, specified mappings between points in an absolute space such as latitude and longitude (φ, λ) and the projected space (x,y). A naive serialization replacing mapstring in an actual proj string might proceed as follows: ′(φ1,λ1),(x1,y1);(φ2,λ2),(x2,y2);...(φN,λN),(xN,yN);′ All φi, λi, xi, and yi are likewise to be replaced by the appropriate real numbers.

To determine how an arbitrary point would be transformed, as suggested above, different interpolation methods could be defined. One such interpolation would involve having first computed the Delaunay triangulation of the known source points. The arbitrary point would be located within one of those triangles and the associated known transformation of that triangle into the output coordinate system would be used to transform the arbitrary point.

We know that proj strings can’t be infinite in practice, for a variety of reasons, but those limitations still enable a tremendous amount of functionality. We were thinking that one could fit the needed code within a single classical .c file used to implement a projection in proj if needed, though if there were other places for useful code, such as for triangulation or spatial indexes or queries, that would be at least as helpful.

All this being said, we are new to proj as contributors, and indeed, are academics not professional software developers, so we wanted to come to the community first before going further. Your thoughts are most appreciated. Further, we want to acknowledge and thank you for all your varied efforts to create, maintain, and extend this key node in the open source geo world!

Best wishes,
Luke


Happy to send you a copy of:

Luke R. Bergmann and David O’Sullivan. 2017. Computing with many spaces: Generalizing projections for the digital geohumanities and GIScience. In Proceedings of GeoHumanities’17: 1st ACM SIGSPATIAL Workshop on Geospatial Humanities, Los Angeles Area, CA, USA, November 7–10, 2017 (GeoHumanities’17), 8 pages.
https://doi.org/10.1145/3149858.3149866

https://dl.acm.org/citation.cfm?id=3149866
Post by Even Rouault
Post by Kristian Evers
Why is that? To me it seems somewhat similar to polynomial mappings
(Horner), in the way that many coefficients are likely to be stored for a
given transformation. Given a clever way of defining the triangles and
their corresponding Helmet parameters I think it should be possible to make
these transformations work. At least as PROJ strings.
You would need to store hundreds of triangles and their parameters, and have
PROJ figure out in which triangle a point it. Certainly doable, but not
directly in the scope of my current work. And avoid the database to grow too
big with very particular transforms. That said, 600 triangles with 6 double
parameters each + the x,y coordinates is just 38.4 KB if you store that as a
compact binary blob. Some balance to know if it is something acceptable in-db
(EPSG style transformations have at most 18 parameters each) or out-db
(grids). The suggestion for the grid approach was to make it more immediatly
usable.
--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
Proj mailing list
http://lists.maptools.org/mailman/listinfo/proj
Kristian Evers
2018-10-29 19:25:11 UTC
Permalink
Luke,

It is an interesting idea you and your colleague are bringing to the table here. I can certainly see why it is intriguing to have data-driven projections readily available in GIS applications. I am not entirely convinced that PROJ is the right tool for the job though. I can easily see PROJ implement various generic mappings consisting of an arbitrary set of parameters [0], where it would be up to the user of PROJ to submit the needed parameters. Similarly to how various projection parameters are in the hands of the user already.

The paper doesn’t go into much details with the math behind the mappings, nor how the mapping between geodetic coordinates and projected coordinates stated in the mapstring should be determined (it seems like there is a chicken and egg type of problem here). For me to really be able make a judgement if this is suitable for PROJ I need some elaboration on those points.

/Kristian

[0] PROJ’s Horner transformation is an example of one such mapping (https://proj4.org/operations/transformations/horner.html?highlight=horner)

On 28 Oct 2018, at 19:54, Luke Bergmann <***@gmail.com<mailto:***@gmail.com>> wrote:
Kristian Evers
2018-10-29 07:36:18 UTC
Permalink
Even,

I wasn't suggesting that this was something you should take on. Your
comment just peaked my interest as I have been pondering this problem
before.

A gridded solution would definitely be the easiest to implement, but not
necessarily the most satisfactory. I for one would always prefer a solution
as close to the original transformation definition as possible.

It would probably also be possible to invent a simple text based file format
for storing the triangles and their related parameters. This would be sort
of equivalent to a grid file that can be linked to in the database. Is the
current setup geared for something like that, or is it strictly reserved for
grid files?

/Kristian

-----Oprindelig meddelelse-----
Fra: Even Rouault <***@spatialys.com>
Sendt: 28. oktober 2018 17:46
Til: ***@lists.maptools.org
Cc: Kristian Evers <***@sdfe.dk>; Andre Joost <andre+***@nurfuerspam.de>
Emne: Re: [Proj] Understanding 2D Helmert
Post by Kristian Evers
Why is that? To me it seems somewhat similar to polynomial mappings
(Horner), in the way that many coefficients are likely to be stored for a
given transformation. Given a clever way of defining the triangles and
their corresponding Helmet parameters I think it should be possible to make
these transformations work. At least as PROJ strings.
You would need to store hundreds of triangles and their parameters, and have
PROJ figure out in which triangle a point it. Certainly doable, but not
directly in the scope of my current work. And avoid the database to grow too
big with very particular transforms. That said, 600 triangles with 6 double
parameters each + the x,y coordinates is just 38.4 KB if you store that as a
compact binary blob. Some balance to know if it is something acceptable in-db
(EPSG style transformations have at most 18 parameters each) or out-db
(grids). The suggestion for the grid approach was to make it more immediatly
usable.
--
Spatialys - Geospatial professional services
http://www.spatialys.com
Andre Joost
2018-10-29 08:44:10 UTC
Permalink
Kristian,

if you register to the website I mentioned earlier, you get those plain
text files, 385 kB uncompressed for both directions. For use in PROJ you
would need an official grant to redistribute them. Looking into the
files, it looks like a 6-parms affine transformation instead of simple
Helmert 2D with 4 parameters.

The main difference to current ntv2-like horizontal grids is that this
method transforms from projected coordinates to projected coordinates,
instead of reprojecting to degrees, transforming, and reprojecting to a
new projected CRS. So you have to code something new anyhow.

Would be a nice idea for GSOC, if someone is interested to dive into it.

Greetings,
Andre Joost
Post by Kristian Evers
Even,
I wasn't suggesting that this was something you should take on. Your
comment just peaked my interest as I have been pondering this problem
before.
A gridded solution would definitely be the easiest to implement, but not
necessarily the most satisfactory. I for one would always prefer a solution
as close to the original transformation definition as possible.
It would probably also be possible to invent a simple text based file format
for storing the triangles and their related parameters. This would be sort
of equivalent to a grid file that can be linked to in the database. Is the
current setup geared for something like that, or is it strictly reserved for
grid files?
/Kristian
-----Oprindelig meddelelse-----
Sendt: 28. oktober 2018 17:46
Emne: Re: [Proj] Understanding 2D Helmert
Post by Kristian Evers
Why is that? To me it seems somewhat similar to polynomial mappings
(Horner), in the way that many coefficients are likely to be stored for a
given transformation. Given a clever way of defining the triangles and
their corresponding Helmet parameters I think it should be possible to make
these transformations work. At least as PROJ strings.
You would need to store hundreds of triangles and their parameters, and have
PROJ figure out in which triangle a point it. Certainly doable, but not
directly in the scope of my current work. And avoid the database to grow too
big with very particular transforms. That said, 600 triangles with 6 double
parameters each + the x,y coordinates is just 38.4 KB if you store that as a
compact binary blob. Some balance to know if it is something acceptable in-db
(EPSG style transformations have at most 18 parameters each) or out-db
(grids). The suggestion for the grid approach was to make it more immediatly
usable.
Even Rouault
2018-10-29 09:58:14 UTC
Permalink
Post by Kristian Evers
It would probably also be possible to invent a simple text based file format
for storing the triangles and their related parameters. This would be sort
of equivalent to a grid file that can be linked to in the database. Is the
current setup geared for something like that, or is it strictly reserved
for grid files?
There are 2 possibilities I can see with the current state of the code. Either
use the table that is dedicated to grid transform and invent a method name for
that method, and add the code that translates this into a PROJ string. This
has the advantage that the method that the CoordinateOperation::gridsNeeded()
can be used and users can know that the transform depends on a external
resource they must possibly grab from somewhere. Or just use the capability of
putting an arbitrary PROJ string to describe a transformation between 2 CRSs
codes (blackbox approach).

Even
--
Spatialys - Geospatial professional services
http://www.spatialys.com
Andre Joost
2018-10-28 17:46:01 UTC
Permalink
Post by Even Rouault
Post by Andre Joost
I have 8 lines within my input file, but only the first 7 get
transformed. The last input line does not end with CR LF.
Would be worth opening a ticket about that.
Done: https://github.com/OSGeo/proj.4/issues/1166

Can you test it on Linux?

Greetings,
Andre Joost
Loading...