# Map projection lat, lon to pixel?

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):

• Do you know what datum your latitude and longitude are in? – Mintx Sep 16 '13 at 21:30
• I think it's spherical earth but I don't know. My projection is spherical mercantor? Where can I find better database? – Gigamegs Sep 16 '13 at 21:45
• Is this question homework? See here for a related question: gis.stackexchange.com/questions/50140/… – Mintx Sep 16 '13 at 21:56
• No. It's not homework.The weekend is yesterday. The link you provided doesn't work because it's for a TILE SERVER. I don't need tiles. I didn't wrote it in my question either. I can't believe it's so difficult but I didn't found many in the internet. Thanks for reading. – Gigamegs Sep 16 '13 at 22:02
• I think it's wgs84:arnulf.us/PLZ – Gigamegs Sep 17 '13 at 8:47

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:

// 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);
}
``````