Coordinate System – How to Map Projection from Lat, Lon to Pixel

Question

I'm trying to find the concave hull of a lat, lon pointset but when I project my points on my map it doesn't look like expected. I've tried a lineare projection and a spherical projection. My data is about 10000-70000 lat lon pairs. How is this possible?

$mapWidth    = 2000;
$mapHeight   = 2000;
$mapLonDelta = $mapLonRight - $mapLonLeft;
$mapLatDelta = $mapLatTop - $mapLatBottom;

$worldMapWidth=(($mapWidth/$mapLonDelta)*360)/(2*M_PI);
$f=min(max(sin($mapLatBottom*(M_PI/180)),-0.9999),0.9999);
$mapOffsetY=$worldMapWidth/2 * log((1+$f)/(1-$f));
$f=min(max(sin($mapLatTop*(M_PI/180)),-0.9999),0.9999);
$mapOffsetTopY=$worldMapWidth/2 * log((1+$f)/(1-$f));
$mapHeightNew=$mapOffsetTopY-$mapOffsetY;
$mapRatioH=$mapHeight/$mapHeightNew;
$newWidth=$mapWidth*($mapHeightNew/$mapHeight);
$mapRatioW=$mapWidth/$newWidth;

$tx = ($lon - $mapLonLeft) * ($newWidth/$mapLonDelta)*$mapRatioW;
$f = sin($lat*M_PI/180);
$ty = ($mapHeightNew-(($worldMapWidth/2 * log((1+$f)/(1-$f)))-$mapOffsetY));

This should look like Germany. What's wrong?:

Update: It seems to work. This is a map of Bonn(Germany):

Answer 1

For files referenced in this answer, see https://gist.github.com/gagern/6636176.

I have some doubts about your input data. Using the data base of postal codes you mentioned, I computed the center of gravity for every polygon to obtain a set of points in Germany. Feeding that data into your algorithm, as posted in the question, I obtained the following image:

Modifying the map dimensions, e.g. $mapWidth = 1000 and $mapHeight = 1500, I was able to see all of Germany. So maybe as a first step would be comparing your set of coordinates to mine. To help with this, I included the postal code for every data point in my data file. I hope this helps you compare them to your data.

Nevertheless, getting parameters right in your code seems a bit tricky. I'd suggest an alternative where the chosen bounding box in geographic coordinates will exactly match the chosen size of the map. So here is the code I'd use:

// This map would show Germany:
$south = deg2rad(47.2);
$north = deg2rad(55.2);
$west = deg2rad(5.8);
$east = deg2rad(15.2);

// This also controls the aspect ratio of the projection
$width = 1000;
$height = 1500;

// Formula for mercator projection y coordinate:
function mercY($lat) { return log(tan($lat/2 + M_PI/4)); }

// Some constants to relate chosen area to screen coordinates
$ymin = mercY($south);
$ymax = mercY($north);
$xFactor = $width/($east - $west);
$yFactor = $height/($ymax - $ymin);

function mapProject($lat, $lon) { // both in radians, use deg2rad if neccessary
  global $xFactor, $yFactor, $west, $ymax;
  $x = $lon;
  $y = mercY($lat);
  $x = ($x - $west)*$xFactor;
  $y = ($ymax - $y)*$yFactor; // y points south
  return array($x, $y);
}