[Proj] Some tests with proj4

Wed Aug 22 23:13:31 PDT 2018

Rethinking, I realize that this is probably not what you want, since the
height coordinate will then be referred to the spherical reference, rather
than the ellipsoidal (hence the 1000 m height in the first example, growing
to 12429 m, when reckoned from the ellipsoid).

Any better suggestions?

Den ons. 22. aug. 2018 kl. 17.26 skrev Thomas Knudsen <
knudsen.thomas at gmail.com>:

>
> 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.osgeo.org/pipermail/proj/attachments/20180823/410aab69/attachment.html>