Mapscript PHP Projection problem

[tbonfort]

Reporter: sven.borgers@kineointeractive.com Date: 2006/08/18 - 19:02

There is a problem with projection in mapscript. When translating screen pixel 
xy coordinates to latlon, it is not possible to have Mapscript re-produce 
projections done by mapserver built in template variables and cs2cs proj4 
comand line app. 

The basic mapserver template variables give a the correct result. This can 
replicated with cs2cs but notwith mapscript php project function. 

The test config has a very simple mapfile with following projection spec and no
projection in the layers:
PROJECTION
  "proj=tmerc"
  "lat_0=49.000000000"
  "lon_0=-2.000000000"
  "k=0.999601"
  "x_0=400000.000"
  "y_0=-100000.000"
  "ellps=airy"
  "towgs84=375,-111,431,0,0,0,0"
END

When using mapserver application to find the latlon for:
x : 533604.343541
y: 180566.045573
The mapserver template MAPLAT, MAPLON variables return:
Lat: -0.076218
Lon: 51.508185
Which, looking these up on MapQuest puts the coords on the right spot.

When using the command line proj4 applications the following results are 
generated:
invproj  -f "%.6f" +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601
+x_0=400000 +y_0=-100000 +ellps=airy +towgs84=375,-111,431,0,0,0,0 <<EOF
> 533604.343541 180566.045573
> EOF
Lat: -0.074626      
Lon: 51.507727
Which is not exactly the same point.

cs2cs command gives following result:
cs2cs  -f "%.6f" +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000
+y_0=-100000 +ellps=airy +towgs84=375,-111,431,0,0,0,0 +to +proj=latlong
+ellps=WGS84 +datum=WGS84 <<EOF
> 533604.343541 180566.045573
> EOF
Lat: -0.076218       
Lon: 51.508185 45.074824
Which is the same as what mapserver returns

Mapscript PHP:
   $projInObj = ms_newProjectionObj( "proj=tmerc,".
  "lat_0=49,".
  "lon_0=-2,".
  "k=0.999601,".
  "x_0=400000,".
  "y_0=-100000,".
  "ellps=airy,".
  "towgs84=375,-111,431,0,0,0,0");

  $projOutObj = ms_newProjectionObj("proj=latlong,  ellps=WGS84,
datum=WGS84");

  $poPoint = ms_newPointObj();
  $poPoint->setXY($_REQUEST["lon"],$_REQUEST["lat"]);
  $poPoint->project($projInObj, $projOutObj);

I get: 
Lat: -0.074619464992752 
Lon: 51.505774438372
Which is close but still about 1km out of where it should land.

[tbonfort]

Author: fwarmerdam Date: 2006/08/18 - 19:39

Sven,

I found the problem!  It is that you can't use the comma delimited format
for the projection string when using towgs84 since there is no way to
know that all the towgs84 arements are part of the same argument.

I changed it to this and it worked:

$projInObj = ms_newProjectionObj( "+proj=tmerc ".
  "+lat_0=49 ".
  "+lon_0=-2 ".
  "+k=0.999601 ".
  "+x_0=400000 ".
  "+y_0=-100000 ".
  "+ellps=airy ".
  "+towgs84=375,-111,431,0,0,0,0");

I would like to strongly encourage folks to use the "+" format for all
MapServer purposes!