[Proj] Some tests with proj4

Wed Aug 22 10:26:56 EST 2018

Hello Jean-Christophe,

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.

proj_geocentric_latitude() works correctly (I believe) *** but only on the
surface of the ellipsoid ***.

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.

I think you will get what you expect if you do as follows:

cct +proj=pipeline +step +proj=cart +R=3396190 +step +proj=cart +inv
+a=3396190 +b=3376200

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.

At least in this case the two latitudes converge as the height grows:

echo 2 49 1e3 0 | cct ...
  2.0000000000 49.3334430922  12429.1901 0.0000

echo 2 49 1e12 0 | cct ...
  2.0000000000   49.0000011347  1000000011372.7307 0.0000

/Thomas

Den ons. 22. aug. 2018 kl. 13.11 skrev Jean-Christophe Malapert <
jcmalapert at gmail.com>:

>
> Hi Kristian
>
> Back to my proj4/gdal problem that I have posted on the proj4 mailing list
>
> I did some tests  and I woud like to share them with you
>
> 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
> echo "-4025114 2971144" | gdaltransform *-s_srs* '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",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]]'
> *-t_srs* 'GEOGCS["Mars 2015 Ocentric",DATUM["D_Mars_2015
> ",SPHEROID["Mars_2015_IAU_Ocentric",3396190,*0*
> ],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],AUTHORITY["IAU2015","49900"]]'
>
> -67.9061078028676 50.1249964005599 0
>
> 2/ Same request but I want to ographic CRS as output
>
> echo "-4025114 2971144" | gdaltransform *-s_srs* '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",*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]]'
> *-t_srs* 'GEOGCS["Mars 2015 Ographic",DATUM["D_Mars_2015
> ",SPHEROID["Mars_2015_IAU_Ographic",3396190,*169.89444722361*
> ],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],AUTHORITY["IAU2015","49901"]]'
>
> -67.9061078028676 50.1249964005599 0
>
> Conclusion : The result is the same and it is not possible !!! - proj4 is
> not intialized correctly in GDAL
>
> 3/ Same request as 2 but I force the init of proj4 with the EXTENSION
> parameter
> echo "-4025114 2971144" | gdaltransform *-s_srs* '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",*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],*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"]*]' *-t_srs* 'GEOGCS["Mars 2015
> Ographic",DATUM["D_Mars_2015 ",SPHEROID["Mars_2015_IAU_Ographic",3396190,
> *169.89444722361*
> ],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],*EXTENSION["PROJ4","+proj=longlat
> +a=3396190 +b=3376200 +towgs84=0,0,0 +no_defs"]*
> ,AUTHORITY["IAU2015","49901"]]'
>
> -67.9061078028676 50.4563269924926 11816.1605372787
>
> conclusion : the latitude changes, much better but is it correct ?
>
> 4/ tranformation ocentric to ographic sans EXTENSION
> echo "0 70" | gdaltransform *-s_srs* 'GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_IAU_IAG_Sphere
> ",3396190,*0*]],PRIMEM["Airy-0",0],UNIT["degree",0.0174532925199433]]'
> *-t_srs* 'GEOGCS["Mars 2015 Ographic",DATUM["D_Mars_2015
> ",SPHEROID["Mars_2015_IAU_Ographic",3396190,*169.89444722361*
> ],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],AUTHORITY["IAU2015","49901"]]'
>
> 0 70 0
>
> Conclusion : same problem than 1
>
> 5/ transformation ocentric to ographic avec EXTENSION
> echo "0 70" | gdaltransform *-s_srs *'GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_IAU_IAG_Sphere
> ",3396190,*0*]],PRIMEM["Airy-0",0],UNIT["degree",0.0174532925199433],*EXTENSION["PROJ4","+proj=longlat
> +a=3396190 +b=3396190 +towgs84=0,0,0 +no_defs"]*]' *-t_srs* 'GEOGCS["Mars
> 2015 Ographic",DATUM["D_Mars_2015
> ",SPHEROID["Mars_2015_IAU_Ographic",3396190,*169.89444722361*
> ],AUTHORITY["IAU","WGCCRE"]],PRIMEM["Airy-0",0,AUTHORITY["IAU","WGCCRE"]],UNIT["Degree",0.017453292519943295],*EXTENSION["PROJ4","+proj=longlat
> +a=3396190 +b=3376200 +towgs84=0,0,0 +no_defs"]*
> ,AUTHORITY["IAU2015","49901"]]'
>
> 0 70.2153180967239 17669.7044292935
>
> conclusion : the latitude changes, much better but is it correct ?
>
> 6/ Even's tests
> $ cs2cs -f "%.8f" +proj=longlat +a=3396190 +b=3396190 +towgs84=0,0,0 +to
> \
>                   +proj=longlat +a=3396190 +b=3376200 +towgs84=0,0,0
> 2 49 0
> 2.00000000      49.33354110 11429.20699246
>
>
> *With PROJ5*
> $ cct +proj=pipeline +step +proj=geoc +inv +a=3396190 +b=3376200
> 2 49 0
>   2.0000000000   49.3346654289        0.0000           inf
>
> 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
>
> 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 *+proj=longlat +a=3396190 +b=3396190 **+no_defs*
> 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.
>
> Is there a way to provide a patch to PROJ4 the time PROJ5 is released with
> the stack gdal/mapserver ?
>
> Thanks a lot,
> Jean-Christophe
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20180822/842987f0/attachment.htm 