[gdal-dev] GDAL 3, PROJ 6.2 exportToProj4() drops Datum silently

[Michael Sumner]

Great discussion, thanks to both of you.  I have a question.

I wonder how we ought to specify a projection that we desire, I work in
contexts where 1.5m or even 1.5 km is not an unreasonable level of
precision, and we get a lot of value out of specifying local projections
(at ocean basin scale) in LAEA, LCC, STERE and others  (OMERC is especially
good for hermisphere-scale migrating birds, very custom and not related to
human activities) - PROJ strings are obviously an easy to specify the
family and centre point of the projection, and we have a lot of workflows
(computer and human) for this. The datum is a bit irrelevant a lot of the
time but we always include it and try to explain it to users.

Are we expected to express a full WKT or WKT2 string for this?   I know we
can write wrappers in some language, but there doesn't seem to be a culture
around crafting custom projections that I can latch onto. Usually I see
people looking for an EPSG code for their region, and as often as not they
end up using UTM because "it's a projection" and no authority has ever done
a survey in their part of the Southern Ocean and published it with EPSG.

Is something like the gdaltransform pipelining syntax the most concise we
can aim for, or do we somehow cache our own local EPSG-like codes?

Thank you.

Cheers, Mike.



On Wed, Oct 30, 2019 at 7:17 AM Even Rouault <even.rouault at spatialys.com>
wrote:

> > The underlying problem is that users may face silent degradation in
> > workflows as they upgrade sf and rgdal to PROJ >= 6 and GDAL >= 3 if they
> > are not prompted to take remedial steps themselves. I realise that users
> > have become too comfortable, but user guidance will be essential to avoid
> > too much noise on SO, twitter, etc.
>
> I do understand your issues. If contributors are willing to spend time
> bringing back more backward compatibility into PROJ exportToProjString()
> (for
> the +datum part I mean. I wouldn't be enthousiastic of attempting to
> recover
> full backward compatibility on the +towgs84 part), we'd probably accept
> reasonable patches . I consciously did a partial implementation on that
> part
> of the code, so that people have an incentive to embrace other options,
> more
> future proof, like WKT or AUTHORITY:CODE, to encode their CRS.
>
> >
> > > If you use directly the WKT representation from this .prj file with
> PROJ
> > > functions that deal with coordinate transformations (or the
> > > OGRSpatialReference object with GDAL functions), things should work
> > > smoothly.
> >
> > WKT is not an option really, though users may wish to look at the WKT2 of
> > an EPSG code. My attempts to feed WKT into importFromWkt() have failed so
> > far, but that could be string handling (I'm treating the whole WKT block
> > as a single text string).
>
> Using WKT should definitely be considered as an option, and investigating
> why
> that fails is in my opinion a more valuable time investment rather than
> trying
> to revive a ghost from the past. WKT is a standard, PROJ strings are not.
>
> > I realise this too, but to permit legacy results to be replicated,
>
> Some legacy behaviour was just plain wrong/inaccurate. Trying to replicate
> it
> in PROJ would be heartbreaking, an extra maintaince burden, and a
> confusion
> for PROJ users because depending on the path they take, they would get
> different results.
>
> Not-so-random example with GDAL 2.4/PROJ 5 with GDA94 to GDA2020
>
> $ gdalsrsinfo EPSG:4283
>
> PROJ.4 : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
>
> OGC WKT :
> GEOGCS["GDA94",
>     DATUM["Geocentric_Datum_of_Australia_1994",
>         SPHEROID["GRS 1980",6378137,298.257222101,
>             AUTHORITY["EPSG","7019"]],
>         TOWGS84[0,0,0,0,0,0,0],
>         AUTHORITY["EPSG","6283"]],
>     PRIMEM["Greenwich",0,
>         AUTHORITY["EPSG","8901"]],
>     UNIT["degree",0.0174532925199433,
>         AUTHORITY["EPSG","9122"]],
>     AUTHORITY["EPSG","4283"]]
>
> $ gdalsrsinfo EPSG:7844
>
> PROJ.4 : +proj=longlat +ellps=GRS80 +no_defs
>
> OGC WKT :
> GEOGCS["GDA2020",
>     DATUM["Geocentric_Datum_of_Australia_2020",
>         SPHEROID["GRS 1980",6378137,298.257222101,
>             AUTHORITY["EPSG","7019"]],
>         AUTHORITY["EPSG","1168"]],
>     PRIMEM["Greenwich",0,
>         AUTHORITY["EPSG","8901"]],
>     UNIT["degree",0.0174532925199433,
>         AUTHORITY["EPSG","9122"]],
>     AUTHORITY["EPSG","7844"]]
>
> Transforming between the above using PROJ string results in a null
> transformation:
>
> $ echo "150 -30" | gdaltransform -s_srs EPSG:4283 -t_srs EPSG:7844 --debug
> on
> OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs
> OGRCT: Target: +proj=longlat +ellps=GRS80 +no_defs
> [...]
> 150 -30 0
>
>
> Same transformation with GDAL + PROJ master:
>
> $ echo "150 -30" | gdaltransform -s_srs EPSG:4283 -t_srs EPSG:7844 --debug
> on
> OGRCT: proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert
> xy_in=deg xy_out=rad step proj=push v_3 step proj=cart ellps=GRS80 step
> proj=helmert x=0.06155 y=-0.01087 z=-0.04019 rx=-0.0394924 ry=-0.0327221
> rz=-0.0328979 s=-0.009994 convention=coordinate_frame step inv proj=cart
> ellps=GRS80 step proj=pop v_3 step proj=unitconvert xy_in=rad xy_out=deg
> step
> proj=axisswap order=2,1 (GDA94 to GDA2020 (1))
> [...]
> 150.000006074915 -29.9999871787747 0
>
> You get a ~1.5 metre shift to the north east as expected.
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev



-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsumner at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20191030/873fc264/attachment.html>