<div dir="ltr"><br>Hello Jean-Christophe,<br><br>No, the result you get from cct is not strictly correct in a 3D world. The problem is that proj=geoc calls the 4D API function proj_geocentric_latitude(), which is copied directly from code in the original pj_fwd / pj_inv code and reflects the 2D/cartographic origin of PROJ.<br><br>proj_geocentric_latitude() works correctly (I believe) *** but only on the surface of the ellipsoid ***.<br><br>In a true 3D world, the planetocentric latitude should converge toward the planetographical as the height grows. But the proj=geoc implementation currently ignores the height entirely.<br><br>I think you will get what you expect if you do as follows:<br><br>cct +proj=pipeline +step +proj=cart +R=3396190 +step +proj=cart +inv +a=3396190 +b=3376200<br><br>i.e. first convert your planetocentric coordinates to 3D cartesian, using a spherical planet model, then convert back to planetographical coordinates using the proper ellipsoidal model.<br><br>At least in this case the two latitudes converge as the height grows:<br><br>echo 2 49 1e3 0 | cct ...<br> 2.0000000000 49.3334430922 12429.1901 0.0000<br><br>echo 2 49 1e12 0 | cct ...<br> 2.0000000000 49.0000011347 1000000011372.7307 0.0000<br><br>/Thomas<br></div><br><div class="gmail_quote"><div dir="ltr">Den ons. 22. aug. 2018 kl. 13.11 skrev Jean-Christophe Malapert <<a href="mailto:jcmalapert@gmail.com">jcmalapert@gmail.com</a>>:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="auto"><div><div class="gmail_quote"><div dir="ltr"><br></div><div dir="ltr"><div>Hi Kristian</div><div><br></div><div>Back to my proj4/gdal problem that I have posted on the proj4 mailing list</div><div><br></div><div>I did some tests and I woud like to share them with you </div><div><br></div><div>1/ I took a pixel coordinate on meters on my map (this is equirecangular projection with a ocentric datum) and I asked to transform in (long,lat) with a ocentric datum </div><div>echo "-4025114 2971144" | gdaltransform <b>-s_srs</b> 'PROJCS["Mars Simple Cylindrical (Ellipse, clon=0)",GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_IAU_IAG_Sphere ",3396190,<b>0</b>]],PRIMEM["Airy-0",0],UNIT["degree",0.0174532925199433]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["standard_parallel_1",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1]]' <b>-t_srs</b> 'GEOGCS["Mars 2015 Ocentric",DATUM["D_Mars_2015 ",SPHEROID["Mars_2015_IAU_Ocentric",3396190,<b>0</b>],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],AUTHORITY["IAU2015","49900"]]'</div><div><br></div><div>-67.9061078028676 50.1249964005599 0<br></div><div><br></div><div>2/ Same request but I want to ographic CRS as output</div><div><br></div><div><div>echo "-4025114 2971144" | gdaltransform <b>-s_srs</b> 'PROJCS["Mars Simple Cylindrical (Ellipse, clon=0)",GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_IAU_IAG_Sphere ",3396190,0]],PRIMEM["Airy-0",<b>0</b>],UNIT["degree",0.0174532925199433]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["standard_parallel_1",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1]]' <b>-t_srs</b> 'GEOGCS["Mars 2015 Ographic",DATUM["D_Mars_2015 ",SPHEROID["Mars_2015_IAU_Ographic",3396190,<b>169.89444722361</b>],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],AUTHORITY["IAU2015","49901"]]'</div><div><br></div></div><div>-67.9061078028676 50.1249964005599 0<br></div><div><br></div><div>Conclusion : The result is the same and it is not possible !!! - proj4 is not intialized correctly in GDAL</div><div><br></div><div>3/ Same request as 2 but I force the init of proj4 with the EXTENSION parameter</div><div><div>echo "-4025114 2971144" | gdaltransform <b>-s_srs</b> 'PROJCS["Mars Simple Cylindrical (Ellipse, clon=0)",GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_IAU_IAG_Sphere ",3396190,0]],PRIMEM["Airy-0",<b>0</b>],UNIT["degree",0.0174532925199433]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["standard_parallel_1",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1],<b>EXTENSION["PROJ4","+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=3396190 +b=3396190 +units=m +towgs84=0,0,0 +no_defs"]</b>]' <b>-t_srs</b> 'GEOGCS["Mars 2015 Ographic",DATUM["D_Mars_2015 ",SPHEROID["Mars_2015_IAU_Ographic",3396190,<b>169.89444722361</b>],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],<b>EXTENSION["PROJ4","+proj=longlat +a=3396190 +b=3376200 +towgs84=0,0,0 +no_defs"]</b>,AUTHORITY["IAU2015","49901"]]'</div><div><br></div></div><div>-67.9061078028676 50.4563269924926 11816.1605372787<br></div><div><br></div><div>conclusion : the latitude changes, much better but is it correct ?</div><div><br></div><div>4/ tranformation ocentric to ographic sans EXTENSION</div><div><div>echo "0 70" | gdaltransform <b>-s_srs</b> 'GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_IAU_IAG_Sphere ",3396190,<b>0</b>]],PRIMEM["Airy-0",0],UNIT["degree",0.0174532925199433]]' <b>-t_srs</b> 'GEOGCS["Mars 2015 Ographic",DATUM["D_Mars_2015 ",SPHEROID["Mars_2015_IAU_Ographic",3396190,<b>169.89444722361</b>],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],AUTHORITY["IAU2015","49901"]]'</div><div><br></div></div><div>0 70 0<br></div><div><br></div><div>Conclusion : same problem than 1</div><div><br></div><div>5/ transformation ocentric to ographic avec EXTENSION</div><div><div>echo "0 70" | gdaltransform <b>-s_srs </b>'GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_IAU_IAG_Sphere ",3396190,<b>0</b>]],PRIMEM["Airy-0",0],UNIT["degree",0.0174532925199433],<b>EXTENSION["PROJ4","+proj=longlat +a=3396190 +b=3396190 +towgs84=0,0,0 +no_defs"]</b>]' <b>-t_srs</b> 'GEOGCS["Mars 2015 Ographic",DATUM["D_Mars_2015 ",SPHEROID["Mars_2015_IAU_Ographic",3396190,<b>169.89444722361</b>],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],<b>EXTENSION["PROJ4","+proj=longlat +a=3396190 +b=3376200 +towgs84=0,0,0 +no_defs"]</b>,AUTHORITY["IAU2015","49901"]]'</div><div><br></div><div>0 70.2153180967239 17669.7044292935<br></div><div><br></div><div>conclusion : the latitude changes, much better but is it correct ?</div><div><br></div><div>6/ Even's tests</div><div>$ <span class="m_-2724348281703035742m_-5612080049211744838gmail-il">cs2cs</span> -f "%.8f" +proj=longlat +a=3396190 +b=3396190 +towgs84=0,0,0 +to \ <br> +proj=longlat +a=3396190 +b=3376200 +towgs84=0,0,0<br>2 49 0<br>2.00000000 49.33354110 11429.20699246<br><br><br><b>With PROJ5</b><br>$ cct +proj=pipeline +step +proj=geoc +inv +a=3396190 +b=3376200<br>2 49 0<br> 2.0000000000 49.3346654289 0.0000 inf<br></div><div><br></div><div>Conclusion : PROJ5 abd PROJ4 give different results. It seems the result of PROJ5 is better because of the third value in the result. It means towgs84=0,0,0 is not enough to get ographic coordinates from ocentric coordinates</div><div><br></div><div>The current problem is that a lot of SIG softwares use PROJ4 (mapserver, gdal, QGIS). For instance QGIS uses Mars2000 CRS with the following proj4 string <b>+proj=longlat +a=3396190 +b=3396190 </b><b>+no_defs</b></div><div>In my case, my data is projected but with an ocentric datum. It means when I think to get a ographic latitude, I get an ocentric latitude.</div><div><br></div><div>Is there a way to provide a patch to PROJ4 the time PROJ5 is released with the stack gdal/mapserver ?</div><div><br></div></div><div>Thanks a lot,</div><div>Jean-Christophe</div></div>
</div></div></div>
_______________________________________________<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/mailman/listinfo/proj</a></blockquote></div>