<div dir="auto">Yes, if datum is not specified on both sides, or if the two datums pass as identical, pj_transform() will not attempt to do any kind of datum shift.<div dir="auto"><br></div><div dir="auto">So do you imply that we have introduced a bug where specifying a non-null grid somehow passes as specifying a datum while specifying the null grid doesn't?</div><div dir="auto"><br></div><div dir="auto">Under all circumstances, I'm not sure how it should be interpreted: we claim to be on NAD83/WGS84 or whatever looks like it at the meter level. We will, however, insist on doing our inverse projection on a spherical, rather than ellipsoidal, reference surface.</div><div dir="auto"><br></div><div dir="auto">In the new API that kind of stuff make sense: We specify what we want done, and libproj does it for us, no matter how stupid it sounds.</div><div dir="auto"><br></div><div dir="auto">But for the old API, libproj is supposed to figure out what to do?</div><div dir="auto"><br></div><div dir="auto">I'm not sure this case makes strict geodetic sense, but somehow libproj used to do the right (for a certain definition of right) thing anyway - and obviously we need to convince it to keep doing that...</div><div dir="auto"><br></div><div dir="auto"><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">Den 1. mar. 2018 8.45 AM skrev "Kristian Evers" <<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>>:<br type="attribution"><blockquote class="quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="word-wrap:break-word;line-break:after-white-space">
I am not sure that is correct, Thomas. See the following output:
<div><br>
</div>
<div>
<div>
<div>$ echo 487147.594520173 4934316.46263998 | ./cs2cs -v +proj=utm +zone=15 +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null</div><div class="quoted-text">
<div># ---- From Coordinate System ----</div>
<div>#Universal Transverse Mercator (UTM)</div>
<div># Cyl, Sph</div>
<div># zone= south</div>
</div><div># +proj=utm +zone=15 +ellps=WGS84</div><div class="quoted-text">
<div># ---- To Coordinate System ----</div>
<div>#Mercator</div>
<div># Cyl, Sph&Ell</div>
<div># lat_ts=</div>
<div># +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0</div>
<div># +units=m +nadgrids=@null</div>
</div><div>-10370728.80 5552839.74 0.00</div>
</div>
<div><br>
</div>
<div>
<div>
<div>
<div>$ echo 487147.594520173 4934316.46263998 | ./cs2cs -v +proj=utm +zone=15 +datum=NAD83 +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m</div><div class="quoted-text">
<div># ---- From Coordinate System ----</div>
<div>#Universal Transverse Mercator (UTM)</div>
<div># Cyl, Sph</div>
<div># zone= south</div>
<div># +proj=utm +zone=15 +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0</div>
<div># ---- To Coordinate System ----</div>
<div>#Mercator</div>
<div># Cyl, Sph&Ell</div>
<div># lat_ts=</div>
<div># +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0</div>
<div># +units=m</div>
</div><div>-10370728.80 5552839.74 0.0</div>
</div>
</div>
</div>
<div><br>
</div>
<div>It seems to be the combination of +datum=NAD83 in the input system and +nadgrids=@null in the output system.</div>
<div>When only one (either one) of them is present the transformation is fine. They are supposed to be treated as being</div>
<div>the same datum, but isn’t. I have an idea where to fix this.</div><font color="#888888">
<div><br>
</div>
<div>/Kristian</div></font><div class="elided-text">
<div><br>
<blockquote type="cite">
<div>On 1 Mar 2018, at 08:15, Thomas Knudsen <<a href="mailto:knudsen.thomas@gmail.com" target="_blank">knudsen.thomas@gmail.com</a>> wrote:</div>
<br class="m_-208034120748076388Apple-interchange-newline">
<div>
<div dir="ltr">
<div>Actually, it seems related to handling of the null grid specificly: When using the conus grid, the 10 km jump into the lithosphere is avoided (although we get the 2D correction we ask for, but really don't want):<br>
<br>
$ echo 487147.594520173 4934316.46263998 0 | cs2cs -f %.16f %FROM% +to +proj=latlong +datun=NAD83 | cs2cs +proj=latlong +datum=NAD83 +to +proj=merc +R=6378137 +nadgrids=conus<br>
<br>
> -10370704.71 5552845.16 -0.00 <br>
<br>
</div>
Note: Height set to negative zero... not sure yet whether this is a hint<br>
<div><br>
<br>
<br>
<br>
<br>
</div>
</div>
<div class="gmail_extra"><br>
<div class="gmail_quote">2018-03-01 7:58 GMT+01:00 Thomas Knudsen <span dir="ltr">
<<a href="mailto:knudsen.thomas@gmail.com" target="_blank">knudsen.thomas@gmail.com</a>></span>:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr">
<div>
<div>
<div>
<div>
<div>Jeff,<br>
<br>
</div>
That's interesting (although annoying!)<br>
<br>
</div>
When I leave out the final nadgrids=@null, I get the correct result from PROJ v5RC6. Can you reproduce that?<br>
</div>
<br>
Note, I'm not implying that this should be considered a solution or workaround, merely trying to narrow in on the<br>
</div>
problem, which seems to be related to datum handling.<span class="m_-208034120748076388HOEnZb"><font color="#888888"><br>
<br>
</font></span></div>
<span class="m_-208034120748076388HOEnZb"><font color="#888888">/Thomas<br>
<div>
<div>
<div><br>
<br>
<div>
<div>
<div><br>
<br>
</div>
</div>
</div>
</div>
</div>
</div>
</font></span></div>
<div class="m_-208034120748076388HOEnZb">
<div class="m_-208034120748076388h5">
<div class="gmail_extra"><br>
<div class="gmail_quote">2018-03-01 0:38 GMT+01:00 Jeff McKenna <span dir="ltr">
<<a href="mailto:jmckenna@gatewaygeomatics.com" target="_blank">jmckenna@gatewaygeomatics.com</a><wbr>></span>:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Update: I never heard back from anyone else about testing PROJ v5 +<br>
MapServer (and GDAL), so in fact I may be the only person to test this<br>
combination so far. Because this is a PROJ mailing list I first give<br>
some more feedback with PROJ utilities, then later will list results<br>
from other software for anyone interested. (in reality I first noticed<br>
the problem in MapServer but then each test in GDAL then PROJ failed as<br>
well, on Windows with PROJ 5.0.0-RC6)<br>
<br>
cs2cs<br>
=====<br>
<br>
I wanted to test PROJ v5 by using a point<br>
(487147.594520173,4934316.4626<wbr>3998) in EPSG:26915 (UTM Zone 15) and<br>
reprojecting that to EPSG:3857 (Google's Mercator).<br>
<br>
with PROJ 4.9.3<br>
---------------<br>
<br>
cs2cs -v +proj=utm +zone=15 +datum=NAD83 +to +proj=merc +a=6378137<br>
+b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m<br>
+nadgrids=@null<br>
# ---- From Coordinate System ----<br>
#Universal Transverse Mercator (UTM)<br>
# Cyl, Sph<br>
# zone= south<br>
# +proj=utm +zone=15 +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0<br>
# ---- To Coordinate System ----<br>
#Mercator<br>
# Cyl, Sph&Ell<br>
# lat_ts=<br>
# +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0<br>
+y_0=0 +k=1.0<br>
# +units=m +nadgrids=@null<br>
487147.594520173 4934316.46263998<br>
<br>
RESULT POINT: -10370728.80 5552839.74 -0.00<br>
<br>
with PROJ 5.0.0-RC6<br>
-------------------<br>
<br>
>cs2cs -v +proj=utm +zone=15 +datum=NAD83 +to +proj=merc +a=6378137<br>
+b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m<br>
+nadgrids=@null<br>
# ---- From Coordinate System ----<br>
#Universal Transverse Mercator (UTM)<br>
# Cyl, Sph<br>
# zone= south<br>
# +proj=utm +zone=15 +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0<br>
# ---- To Coordinate System ----<br>
#Mercator<br>
# Cyl, Sph&Ell<br>
# lat_ts=<br>
# +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0<br>
+y_0=0 +k=1.0<br>
# +units=m +nadgrids=@null<br>
487147.594520173 4934316.46263998<br>
<br>
RESULT POINT: -10370728.80 5522830.15 -10484.03<br>
<br>
Discussion<br>
----------<br>
<br>
PROJ 4.9.3 returns the correct value (I've confirmed this with many<br>
other tools including QGIS). Can someone explain why those 2 results<br>
are different? (I have reproduced this problem on Windows and also<br>
Ubuntu 14.94)<br>
<br>
GDAL<br>
====<br>
<br>
I wanted to test PROJ v5 by using that same test point (as a shapefile)<br>
in EPSG:26915 (UTM Zone 15) and reprojecting that to EPSG:3857 (Google's<br>
Mercator), using GDAL's ogr2ogr command. (data used:<br>
<a href="https://www.gatewaygeomatics.com/dl/parcels-proj-v5-test.zip" rel="noreferrer" target="_blank">https://www.gatewaygeomatics.c<wbr>om/dl/parcels-proj-v5-test.zip</a> ) Below I<br>
have used the ogr2ogr command (included in the "commands.txt" file in<br>
that data download) and then ogrinfo on the reprojected file to output:<br>
<br>
with PROJ 4.9.3<br>
---------------<br>
<br>
POINT (-10370728.7959208 5552839.74177463)<br>
<br>
with PROJ 5.0.0-RC6<br>
-------------------<br>
<br>
POINT (-10370728.7959208 5522830.1464459)<br>
<br>
Discussion<br>
----------<br>
<br>
PROJ 4.9.3 compiled into GDAL 2.2.3 returns the correct point value, but<br>
returns a different point value with PROJ 5.0.0-RC6<br>
<br>
MapServer<br>
=========<br>
<br>
Same problems as GDAL (reprojected layer has different extents with PROJ<br>
5.0.0-RC6, causing problems to display the reprojected layer in a map<br>
image).<br>
<br>
<br>
If you made it this far in my long email, thanks for listening (I am<br>
totally stumped on this one, again ha).<br>
<span class="m_-208034120748076388m_-1307082090562542093HOEnZb"><font color="#888888"><br>
-jeff<br>
</font></span>
<div class="m_-208034120748076388m_-1307082090562542093HOEnZb">
<div class="m_-208034120748076388m_-1307082090562542093h5"><br>
<br>
<br>
<br>
<br>
On 2018-02-25 3:45 PM, Kristian Evers wrote:<br>
> Thanks for testing, Jeff. It is very appreciated. I hope you sort of your problems,<br>
> PROJ-related or not.<br>
><br>
> /Kristian<br>
><br>
>> On 25 Feb 2018, at 00:29, Jeff McKenna <<a href="mailto:jmckenna@gatewaygeomatics.com" target="_blank">jmckenna@gatewaygeomatics.com</a><wbr>> wrote:<br>
>><br>
>> Update: it seems many things were happening, causing me to look in wrong<br>
>> places. But in the end MapServer was in fact reprojecting the points;<br>
>> what caused me total confusion was that I was testing with several<br>
>> MapServer-related applications that are part of MS4W (Mapbender,<br>
>> GeoMOOSE) and one application had a totally blank opening map image with<br>
>> PROJ v5. In fact for some reason the old hardcoded map extents<br>
>> (displayed in EPSG:3857) no longer display the data (you have to now<br>
>> zoom out a bit to see the data, reprojected with PROJ v5). I'll contact<br>
>> the GeoMOOSE team to let them dive into this (still, it could be my<br>
>> Windows environment, as no one else seems to be testing MapServer + PROJ<br>
>> v5 yet).<br>
>><br>
>> Very sorry for all this noise!<br>
>><br>
>> Short story: all is fine with RC5 here.<br>
>><br>
>> Thank you again Kristian for your passion driving this.<br>
>><br>
>> -jeff<br>
>><br>
>><br>
>><br>
>><br>
>><br>
>> On 2018-02-23 5:33 PM, Jeff McKenna wrote:<br>
>>> Hi Kristian,<br>
>>><br>
>>> No problems here on Windows in the MS4W environment, but I'm running<br>
>>> into odd errors when I have MapServer call PROJ 5<br>
>>> (<a href="https://lists.osgeo.org/pipermail/mapserver-dev/2018-February/015347.html" rel="noreferrer" target="_blank">https://lists.osgeo.org/piper<wbr>mail/mapserver-dev/2018-Februa<wbr>ry/015347.html</a>).<br>
>>> Obviously an issue downstream not here, but I thought I'd at least<br>
>>> mention it here that I am testing RC5 as well :)<br>
>>><br>
>>> Thanks for pushing this release along.<br>
>>><br>
>>> -jeff<br>
>>><br>
>>><br>
>>><br>
>>><br>
>>> On 2018-02-23 3:19 PM, Kristian Evers wrote:<br>
>>>> All,<br>
>>>><br>
>>>> I have prepared PROJ 5.0.0RC5. Downloads are here:<br>
>>>><br>
>>>><br>
>>>> <a href="http://download.osgeo.org/proj/proj-5.0.0RC5.tar.gz" rel="noreferrer" target="_blank">
http://download.osgeo.org/proj<wbr>/proj-5.0.0RC5.tar.gz</a> (<a href="http://download.osgeo.org/proj/proj-5.0.0RC5.tar.gz.md5" rel="noreferrer" target="_blank">http://download.osgeo.org/pro<wbr>j/proj-5.0.0RC5.tar.gz.md5</a>)<br>
>>>> <a href="http://download.osgeo.org/proj/proj-5.0.0RC5.zip" rel="noreferrer" target="_blank">
http://download.osgeo.org/proj<wbr>/proj-5.0.0RC5.zip</a> (<a href="http://download.osgeo.org/proj/proj-5.0.0RC5.zip.md5" rel="noreferrer" target="_blank">http://download.osgeo.org/pro<wbr>j/proj-5.0.0RC5.zip.md5</a>)<br>
>>>><br>
>>>><br>
>>>> Since RC4 problems with vertical grid installation, the Horner<br>
>>>> transformation and faulty polygon area calculation<br>
>>>> in geodesic library has been fixed. Hopefully this will be the final<br>
>>>> release candidate. I will call for a vote for<br>
>>>> promotion to final release on the 27th. Unless serious issues are<br>
>>>> discovered in the mean time.<br>
>>>><br>
>>>> /Kristian<br>
>>>><br>
>>>><br>
______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org" target="_blank">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">http://lists.maptools.org/mail<wbr>man/listinfo/proj</a><br>
</div>
</div>
</blockquote>
</div>
<br>
</div>
</div>
</div>
</blockquote>
</div>
<br>
</div>
______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org" target="_blank">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" target="_blank">http://lists.maptools.org/<wbr>mailman/listinfo/proj</a></div>
</blockquote>
</div>
<br>
</div></div>
</div>
<br>______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">http://lists.maptools.org/<wbr>mailman/listinfo/proj</a><br></blockquote></div><br></div>