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