2

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?: Map

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

  • 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
3

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:

Image of postal codes using original code

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.