[mapserver-users] please check my formulae

[Jan Mantkowski]

Hi!

At
http://www.npm.ac.uk/rsdas/doc/description.html
there are s vew very handy formulas to do lat/lon to xy math:

How to find a latitude/longitude point within an image

The images we produce are in Mercator projection and there are a set of
fairly simple formulae for converting latitude/longitude (lat/lon) positions
to x/y pixel co-ordinates.

rows - number of rows in image
cols - number of columns in image
minlon/maxlon - minimum/maximum longitude of image (in decimal degrees)
minlat/maxlat - minimum/maximum latitude of image (in decimal degrees)
Note that these are the coordinates of the CENTRE of the corner pixels.
DEGTORAD - conversion from degrees to radians (PI/180.0)
ln - natural log

Calculate Y position.

   Ymin = ln (tan (DEGTORAD * (45.0 + (minlat / 2.0))))
   Ymax = ln (tan (DEGTORAD * (45.0 + (maxlat / 2.0))))
   Yint = ln (tan (DEGTORAD * (45.0 + (lat / 2.0))))
   FractY = (Yint - YMin) / (YMax - YMin)
   y = (rows - 1) * (1.0 - FractY)

Calculate X position.

   FractX = (lon - minlon) / (maxlon - minlon)
   x = (cols - 1) * FractX


I made this php script out of it:

$Ymin = log (tan (deg2rad (45 + ($minlat / 2))));
$Ymax = log (tan (deg2rad (45 + ($maxlat / 2))));
$Yint = log (tan (deg2rad (45 + ($lat / 2))));
$y = ($width - 1)*(1 -(($Yint - $YMin)/($YMax - $YMin)));

$x = round(($height - 1) * (($lon - $minlon) / ($maxlon - $minlon)));

Unfortunately something seems to be wrong. I never get usefull results.
Please can you tell me if the formula is right at all?


Thanks,

Jan