Transformation Problem with Merchich Sahara Projections

[ysid]

When using cs2cs to transform into a Merchich Sahara projection (EPSG 26194 or 26195) in OSGeo4W v1.9.1, results can be incorrect. Results are correct when no GDAL/PROJ/OSGeo variables are specified in the environment, but are incorrect when specified. To be clear, I will present an example:

Input coordinate:

• Projection: Geographic WGS84 (EPSG 4326)
• Proj4: +proj=longlat +datum=WGS84 +pm=greenwich +ellps=WGS84
• Coordinate: -16.6666666666667 18.6666666666667

Output coordinate:

• Projection: Merchich Sahara South (EPSG 26195)
• Proj4: +proj=lcc +units=m +towgs84=31,146,47,0,0,0,0 +pm=greenwich +a=6378249.2 +b=6356515 +lon_0=-5.4 +x_0=1500000 +y_0=400000 +lat_0=22.5 +lat_1=22.5 +k_0=0.999616437
• Coordinate: 310159.282196 20443.0660249

Note: Instead of the "+b" value, an alternative "+rf" variable can be specified, and will give results that are virtually identical (rf value is 293.466021293627).

The abovementioned coordinate is the correct value, as verified by the website:

​http://cs2cs.mygeodata.eu/

Similar results are also obtained in PCI and ArcGIS.

When running the above transformation using cs2cs, the following command is used:

cs2cs.exe -f "%.03f%" +proj=longlat +datum=WGS84 +pm=greenwich +ellps=WGS84 +to +proj=lcc +units=m +towgs84=31,146,47,0,0,0,0 +pm=greenwich +a=6378249.2 +rf=293.466021293627 +lon_0=-5.4 +x_0=1500000 +y_0=400000 +lat_0=22.5 +lat_1=22.5 +k_0=0.999616437 input.txt > output.txt

The result depends on the environment setting. When no environment setting specific to PROJ4 or GDAL are set, the results are correct. However, when environment variables are set, the result is very wrong.

• Without environment settings: 310159.282 20443.066 -80.008
• With environment settings: 296432.545 39220.367 -80.008

Interestingly, the ellipsoidal height is the same in both cases.

Environment settings used are typical for an OSGeo4W/GDAL installation, as folllows:

GDAL_DATA: <OSGeo4W_install_dir>\share\gdal
GDAL_DRIVER_PATH: <OSGeo4W_install_dir>\bin\gdalplugins\1.9
GEOTIFF_CSV: <OSGeo4W_install_dir>\share\gdal
OSGEO4W_ROOT: <OSGeo4W_install_dir>
PROJ_LIB: <OSGeo4W_install_dir>\share\proj

Strange results are also seen when using gdalwarp with the same from and to projections and environment variables.

Something in the OSGeo4W installation is causing the coordinates to be transformed incorrectly.

[jlacroix]

This problem comes from a default PROJ4 setting. By default, PROJ4 tries to define the +lat_1 and +lat_2 parameters to work with a U.S. map for all lcc and aea based projection. This is defined in <OSGeo4W_install_dir>\share\proj\proj_def.dat

When not using the environment variables, cs2cs simply didn't find this default configuration file and the good result was produced.

To prevent the use of this default configuration file, PROJ4 has the +no_defs parameters. If you use this parameter, like it is in the epsg file and the cs2cs website, you should get the good result.

cs2cs.exe -f "%.03f%" +proj=longlat +datum=WGS84 +pm=greenwich +ellps=WGS84 +to +proj=lcc +units=m +towgs84=31,146,47,0,0,0,0 +pm=greenwich +a=6378249.2 +rf=293.466021293627 +lon_0=-5.4 +x_0=1500000 +y_0=400000 +lat_0=22.5 +lat_1=22.5 +k_0=0.999616437 +no_defs input.txt

310159.282 20443.066 -80.008

[ysid]

Thanks for the info. I was able to solve the problem by deleting the defaults for <lcc> from the file. May I respectfully suggest that projection defaults are a rather obscure and outdated feature of PROJ4. They probably cause more problems than they solve.

[etourigny]

you should probably mention this in the proj4 list or trac - this here is for osgeo4w which is "just" a packaging of proj4.

[warmerdam]

Yusuf,

Point noted, but I am hesitant to change the behavior at this point.

I believe autogenerated dictionaries like the epsg one all explicitly set these parameter values so it does not affect them.