Add geographic offset transformation method. #1138

rouault · 2018-09-29T08:06:20Z

The Geographic offsets transformation adds an offset to the geographic longitude,
latitude coordinates, and an offset to the ellipsoidal height.
This method is normally only used when low accuracy is tolerated. It is documented
as coordinate operation method code 9619 (for geographic 2D) and 9660 (for
geographic 3D) in the EPSG dataset.

It can also be used to implement the method Geographic2D with Height Offsets
(code 9618) by noting that the input vertical component is a gravity-related
height and the output vertical component is the ellispoid height (dh being
the geoid undulation).

It can also be used to implement the method Vertical offset (code 9616)

It is used for example to transform:

from the old Greek geographic 2D CRS to the newer GGRS87 CRS
from Tokyo + JSLD69 height to WGS 84
from Baltic 1977 height to Black Sea height

It is also useful to document the implicit zero-offset transformation
we do in pipelines such as

+proj=pipeline +step +inv +proj=longlat +ellps=A
               +step +proj=longlat +ellps=B

that can be explicited as

+proj=pipeline +step +inv +proj=longlat +ellps=A
               +step +proj=geogoffset [+dlon=0 +dlat=0 +dh=0]
               +step +proj=longlat +ellps=B

kbevers

This seems like a no-brainer to include. I have a few inlined comments but apart from those I think it looks good.

Regarding your concluding explicit examples, it implies that you move from coordinates referenced on one ellipsoid to another. It actually doesn't. So your final pipeline ends up as a chain of three no-ops. Surely it must be fine to just use one of the three steps?

In case it is useful to someone, this is how to move coordinates from one ellipsoid to another:

+proj=pipeline +step +cart +ellps=A
			   +step +cart +ellps=B +inv

kbevers · 2018-09-29T13:59:01Z

docs/source/operations/transformations/geogoffset.rst

+| **output type**     | Geodetic coordinates (horizontal), meters (vertical)     |
+---------------------+----------------------------------------------------------+
+
+This method is normally only used when low accuracy is tolerated. It is documented


A few points regarding this section of the documentation:

Add a proper reference to EPSG guidance note 7-2 (that's where this is from, right?)

Don't mention the specific method numbers. They don't add much value to me as a reader and they may get changed upstream which we are not tracking systematically.

2\. Don't mention the specific method numbers. They don't add much value to me as a reader and they may get changed upstream which we are not tracking systematically.

I hope they won't !!!! It is quite critical for my current work that the remain stable.

1. Add a proper reference to EPSG guidance note 7-2 (that's where this is from, right?)

Yes, done

I hope they won't !!!! It is quite critical for my current work that the remain stable.

Do they offer any guarantees for that? I realize that it's not likely to happen, I'm just wondering if there's something we can do to prevent unpleasant surprises.

Do they offer any guarantees for that? I realize that it's not likely to happen, I'm just wondering if there's something we can do to prevent unpleasant surprises.

My understanding is that they won't reuse numbers. Those are primary keys in a table. They might deprecate entries in favor of new ones. But normally an abandonned method number will not be reused for something different. My WKT -> PROJ converter checks for either the method name (can be a bit fragile too in case of misspellings...) or method code

kbevers · 2018-09-29T14:03:28Z

docs/source/operations/transformations/geogoffset.rst

+latitude coordinates, and an offset to the ellipsoidal height.
+
+---------------------+----------------------------------------------------------+
+| **Alias**           | geogoffset                                               |


It doesn't matter much but I am not a big fan of this name. Having a hard time to come up with something better though, so this may have to do. Maybe just "offset" or perhaps "goffset"?

There are a lot of potential offsetting methods, some much more complicated. Better be explicit

kbevers · 2018-09-29T14:09:16Z

src/PJ_geogoffset.c

+    PJ_COORD point = {{0,0,0,0}};
+
+    /* offset coordinate */
+    point.lpz.lam = xyz.x - Q->dlon;


If you instead use the 4D function as the main function you can make your intentions a bit clearer by using lpz from the PJ_COORD:

static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { obs.lpz.lam -= Q->dlon; obs.lpz.phi -= Q->phi; obs.lpz.z -= Q->dh; return obs; }

instead of having it look like you subtract geodetic coordinates from projected coordinates as you are forced to in the 3D function.

rouault · 2018-09-29T15:03:37Z

Regarding your concluding explicit examples, it implies that you move from coordinates referenced on one ellipsoid to another. It actually doesn't.

Well it moves from one ellipsoid to another, by doing nothing. Which is a valid option when you don't really know how datum relates to each other (going through geog -> cartesian -> geog isn't necessarily the right thing to do). This is just a matter of documenting that you don't do anything, but yes this is a no-op.

busstoptaktik · 2018-09-29T15:14:02Z

Well it moves from one ellipsoid to another, by doing nothing.

Not really ready to comment on a Saturday afternoon, but IMHO, to move from one ellipsoid to another you should first transform to cartesian 3D (pj_cart), using the parameters of the original ellipsoid, then back using the ellipsoidal parameters of the new ellipsoid.

Regarding naming of the method: I would suggest generalizing it a bit, perhaps under the name "stretch", giving optional parameters +scale=sx,sy,sz,st an +offset=sx,sy,sz,st, defaulting all scales to 1 and all offsets to 0.

Actually, I think it may fit very well into the framework of pj_unitconvert.c

rouault · 2018-09-29T15:39:24Z

I would suggest generalizing it a bit

I was considering this at first, but EPSG Guidance 7-2 distinguishes other linear/affine transformation methods such as "Similarity transformation (9621) ( see https://www.epsg-registry.org/export.htm?wkt=urn:ogc:def:coordinateOperation:EPSG::1035 for an example ) or "Affine parametric transformation" (9624) ( https://www.epsg-registry.org/export.htm?wkt=urn:ogc:def:coordinateOperation:EPSG::10087 ), which involve rotation/skewing and take different input parameters. So it might be difficult for users to express EPSG-like transforms as PROJ strings if we have a single generic way

it may fit very well into the framework of pj_unitconvert.c

Wouldn't be obvious at all to users that unitconvert can do affine transformations as well

kbevers · 2018-09-29T15:55:25Z

Wouldn't be obvious at all to users that unitconvert can do affine transformations as well

Yes, I agree. I made a ticket for an affine transformation a long time ago. I did think about suggesting going down that route as a generic option for this operation but I think it is overkill.

@rouault do you need an affine transformation sometime in the near future for your current work? If so, I can add it.

busstoptaktik · 2018-09-29T15:59:26Z

Wouldn't be obvious at all to users that unitconvert can do affine transformations as well

You should think of it the other way round: What unitconvert actually does is affine transformations.

What I suggest is to make this part of a generic framework with selectable defaults (it's easy to make aliases for already existing operations and only slighly less easy to make different defaults).

rouault · 2018-09-29T18:12:18Z

In the general case, to model any affine transform, we transform from (Xsrc, Ysrc, Zsrc, Tsrc) to (Xdst, Ydst, Zdst) with 14 parameters: tx, ty, tz, sji (i=x,y,z j=x,y,z), tt, scale_t (I don't think that mixing T with other dimensions makes much sense)

Xdst = tx + Xsrc * sxx + Ysrc * sxy + Zsrc * sxz
Ydst = ty + Xsrc * syx + Ysrc * syy + Zsrc * syz
Zdst = tz + Xsrc * szx + Ysrc * szy + Zsrc * szz
Tdst = tt + Tsrc * st

So should we go to a 'affine' method like that ? Do 'tx', 'tx', 'tz', 'sxx', ... 'szz', 'tt', 'st' sound like good parameter names ? And have all offsets default to 0, st = 1, diagonal coefficients = 1, non-diagonal coefficients = 0, so that by default this is the identify transform.

For my use case of Geographic offsets, the user would only specify tx, ty, tz

(Actually (non time-dependent) Helmert is a particular case of this...)

kbevers · 2018-09-29T18:58:19Z

So should we go to a 'affine' method like that ?

It's more generic. I prefer that. We may have to add some unit conversion to the parameters though. At least your case in this PR uses arc seconds.

In Do 'tx', 'tx', 'tz', 'sxx', ... 'szz', 'tt', 'st' sound like good parameter names ?

In #535 I suggested using

+proj=affine +A=... +B=... +C=... +D=... +E=... +F=... +G=... +H=... +I=...
                    +xoff=... +yoff=... +zoff=...

based on examples from Shapely that I thought looked nice. Your approach of simply indexing the matrix may be easier to mentally grasp when converting a matrix into parameters. Instead of x,y,z indicies 1,2,3 could also be used: `s11, s12, s13, s21, ..., s33'.

For my use case of Geographic offsets, the user would only specify tx, ty, tz

(Actually (non time-dependent) Helmert is a particular case of this...)

Yes, but in it's current state it won't work for your application since proj=helmert only accepts cartesian coordinates.

rouault · 2018-09-29T19:41:22Z

Yes, but in it's current state it won't work for your application since proj=helmert only accepts cartesian coordinates.

Actually I'll turn the methd to accept PJ_IO_UNITS_WHATEVER as input and output. For geographic case only offsetting makes sense. Scaling only makes sense for cartesian/projected coordinate space.

rouault · 2018-09-29T19:49:19Z

Wondering if I shouldn't accept dlon and dlat parameters as well, that are equivalent to xoff and yoff, except they are expressed in arc-seconds (like in the EPSG dataset), whereas xoff and yoff would be in the natural unit of the coordinate system, ie radians in the geographic case.

The Geographic offsets transformation adds an offset to the geographic longitude, latitude coordinates, and an offset to the ellipsoidal height. This method is normally only used when low accuracy is tolerated. It is documented as coordinate operation method code 9619 (for geographic 2D) and 9660 (for geographic 3D) in the EPSG dataset. It can also be used to implement the method Geographic2D with Height Offsets (code 9618) by noting that the input vertical component is a gravity-related height and the output vertical component is the ellispoid height (dh being the geoid undulation). It can also be used to implement the method Vertical offset (code 9616) It is used for example to transform: - from the old Greek geographic 2D CRS to the newer GGRS87 CRS - from Tokyo + JSLD69 height to WGS 84 - from Baltic 1977 height to Black Sea height It is also useful to document the implicit zero-offset transformation we do in pipelines such as +proj=pipeline +step +inv +proj=longlat +ellps=A +step +proj=longlat +ellps=B that can be explicited as +proj=pipeline +step +inv +proj=longlat +ellps=A +step +proj=geogoffset [+dlon=0 +dlat=0 +dh=0] +step +proj=longlat +ellps=B

rouault · 2018-10-01T10:39:32Z

OK, I've added a more general affine method, and make geogoffset a particular case of it, mostly for the conveniency of being able to add the lon and lat offsets as arc-seconds instead of radians

busstoptaktik · 2018-10-01T14:27:52Z

docs/source/operations/transformations/geogoffset.rst

+
+Conversion from Baltic 1977 height to Black Sea height::
+
+    proj=geogoffset dh=0.4


An awesome example of geodetic humour :-)

Hum, I'm afraid I didn't grasp the humour part of it... I just borrowed this randomly from EPSG.

Actually, I now see that this may actually have practical applications:

https://belarusdigest.com/story/poland-belarus-and-ukraine-to-revive-the-baltic-black-sea-water-route/

But certainly also a humourous twist...

The joke went over my head too... But yes, this is certainly useful, for instance when converting between local height systems and the upcoming IHRS which, as far as I understand, in most cases will be a simply offset like displayed here.

kbevers · 2018-10-01T19:21:12Z

Looks good to me in it's current state. Any chance we'll get some documentation for the affine transformation as well? :)

…ar case of it (fixes OSGeo#535)

rouault · 2018-10-01T19:32:46Z

Any chance we'll get some documentation for the affine transformation as well? :)

Ah I wrote it but forgot to git add it. Now pushed

rouault force-pushed the geogoffset branch 4 times, most recently from 427792b to 6fb8955 Compare September 29, 2018 12:05

kbevers reviewed Sep 29, 2018

View reviewed changes

rouault force-pushed the geogoffset branch from 6fb8955 to 7c9a377 Compare September 29, 2018 15:23

rouault force-pushed the geogoffset branch 2 times, most recently from 88ef03e to 6e39f95 Compare October 1, 2018 10:33

busstoptaktik reviewed Oct 1, 2018

View reviewed changes

Add a affine transformation method, and make geogoffset as a particul…

fe29f8a

…ar case of it (fixes OSGeo#535)

rouault force-pushed the geogoffset branch from 6e39f95 to fe29f8a Compare October 1, 2018 19:32

kbevers merged commit cc33c13 into OSGeo:master Oct 1, 2018

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add geographic offset transformation method. #1138

Add geographic offset transformation method. #1138

rouault commented Sep 29, 2018 •

edited

Loading

kbevers left a comment

kbevers Sep 29, 2018

rouault Sep 29, 2018

rouault Sep 29, 2018

kbevers Sep 29, 2018

rouault Sep 29, 2018

kbevers Sep 29, 2018

rouault Sep 29, 2018

kbevers Sep 29, 2018

rouault Sep 29, 2018

rouault commented Sep 29, 2018

busstoptaktik commented Sep 29, 2018

rouault commented Sep 29, 2018

kbevers commented Sep 29, 2018 •

edited

Loading

busstoptaktik commented Sep 29, 2018

rouault commented Sep 29, 2018

kbevers commented Sep 29, 2018

rouault commented Sep 29, 2018

rouault commented Sep 29, 2018

rouault commented Oct 1, 2018

busstoptaktik Oct 1, 2018

rouault Oct 1, 2018

busstoptaktik Oct 1, 2018

kbevers Oct 1, 2018

kbevers commented Oct 1, 2018

rouault commented Oct 1, 2018


		Conversion from Baltic 1977 height to Black Sea height::

		proj=geogoffset dh=0.4

Add geographic offset transformation method. #1138

Add geographic offset transformation method. #1138

Conversation

rouault commented Sep 29, 2018 • edited Loading

kbevers left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rouault commented Sep 29, 2018

busstoptaktik commented Sep 29, 2018

rouault commented Sep 29, 2018

kbevers commented Sep 29, 2018 • edited Loading

busstoptaktik commented Sep 29, 2018

rouault commented Sep 29, 2018

kbevers commented Sep 29, 2018

rouault commented Sep 29, 2018

rouault commented Sep 29, 2018

rouault commented Oct 1, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kbevers commented Oct 1, 2018

rouault commented Oct 1, 2018

rouault commented Sep 29, 2018 •

edited

Loading

kbevers commented Sep 29, 2018 •

edited

Loading