[Proj] Some tests with proj4

Wed Aug 22 04:09:36 PDT 2018

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20180822/72737aef/attachment.html>