Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Convert latitude/longitude point to a pixels (x,y) on mercator projection

I'm trying to convert a lat/long point into a 2d point so that I can display it on an image of the world-which is a mercator projection.

I've seen various ways of doing this and a few questions on stack overflow-I've tried out the different code snippets and although I get the correct longitude to pixel, the latitude is always off-seems to be getting more reasonable though.

I need the formula to take into account the image size, width etc.

I've tried this piece of code:

double minLat = -85.05112878; double minLong = -180; double maxLat = 85.05112878; double maxLong = 180;  // Map image size (in points) double mapHeight = 768.0; double mapWidth = 991.0;  // Determine the map scale (points per degree) double xScale = mapWidth/ (maxLong - minLong); double yScale = mapHeight / (maxLat - minLat);  // position of map image for point double x = (lon - minLong) * xScale; double y = - (lat + minLat) * yScale;  System.out.println("final coords: " + x + " " + y); 

The latitude seems to be off by about 30px in the example I'm trying. Any help or advice?

Update

Based on this question:Lat/lon to xy

I've tried to use the code provided but I'm still having some problems with latitude conversion, longitude is fine.

int mapWidth = 991; int mapHeight = 768;  double mapLonLeft = -180; double mapLonRight = 180; double mapLonDelta = mapLonRight - mapLonLeft;  double mapLatBottom = -85.05112878; double mapLatBottomDegree = mapLatBottom * Math.PI / 180; double worldMapWidth = ((mapWidth / mapLonDelta) * 360) / (2 * Math.PI); double mapOffsetY = (worldMapWidth / 2 * Math.log((1 + Math.sin(mapLatBottomDegree)) / (1 - Math.sin(mapLatBottomDegree))));  double x = (lon - mapLonLeft) * (mapWidth / mapLonDelta); double y = 0.1; if (lat < 0) {     lat = lat * Math.PI / 180;     y = mapHeight - ((worldMapWidth / 2 * Math.log((1 + Math.sin(lat)) / (1 - Math.sin(lat)))) - mapOffsetY); } else if (lat > 0) {     lat = lat * Math.PI / 180;     lat = lat * -1;     y = mapHeight - ((worldMapWidth / 2 * Math.log((1 + Math.sin(lat)) / (1 - Math.sin(lat)))) - mapOffsetY);     System.out.println("y before minus: " + y);     y = mapHeight - y; } else {     y = mapHeight / 2; } System.out.println(x); System.out.println(y); 

When using the original code if the latitude value is positive it returned a negative point, so I modified it slightly and tested with the extreme latitudes-which should be point 0 and point 766, it works fine. However when I try a different latitude value ex: 58.07 (just north of the UK) it displays as north of Spain.

like image 414
drunkmonkey Avatar asked Jan 15 '13 01:01

drunkmonkey


People also ask

How do you calculate Mercator projection?

The formula for Mercator's projection is T(ϕ, θ)=(θ, ln(|sec(ϕ) + tan(ϕ)|)). Of course, there are a huge number of map projections.

Is latitude considered X or Y?

In this case, longitude values are considered the x-coordinate, while latitude values are the y-coordinate. A geographic coordinate system has the following components: Angular units: The unit of measure on the spherical reference system.

How do you change XY coordinates to latitude and longitude?

Calculate latitude and longitude using the formula: latitude = asin (z/R) and longitude = atan2 (y,x). In this formula, we have the values of x, y, z and R from step 2. Asin is arc sin, which is a mathematical function, and atan2 is a variation of the arc tangent function. The symbol * stands for multiplication.


2 Answers

The Mercator map projection is a special limiting case of the Lambert Conic Conformal map projection with the equator as the single standard parallel. All other parallels of latitude are straight lines and the meridians are also straight lines at right angles to the equator, equally spaced. It is the basis for the transverse and oblique forms of the projection. It is little used for land mapping purposes but is in almost universal use for navigation charts. As well as being conformal, it has the particular property that straight lines drawn on it are lines of constant bearing. Thus navigators may derive their course from the angle the straight course line makes with the meridians. [1.]

Mercator projection

The formulas to derive projected Easting and Northing coordinates from spherical latitude φ and longitude λ are:

E = FE + R (λ – λₒ) N = FN + R ln[tan(π/4 + φ/2)]    

where λO is the longitude of natural origin and FE and FN are false easting and false northing. In spherical Mercator those values are actually not used, so you can simplify the formula to

derivation of the mercator projection (wikipedia)

Pseudo code example, so this can be adapted to every programming language.

latitude    = 41.145556; // (φ) longitude   = -73.995;   // (λ)  mapWidth    = 200; mapHeight   = 100;  // get x value x = (longitude+180)*(mapWidth/360)  // convert from degrees to radians latRad = latitude*PI/180;  // get y value mercN = ln(tan((PI/4)+(latRad/2))); y     = (mapHeight/2)-(mapWidth*mercN/(2*PI)); 

Sources:

  1. OGP Geomatics Committee, Guidance Note Number 7, part 2: Coordinate Conversions and Transformation
  2. Derivation of the Mercator projection
  3. National Atlas: Map Projections
  4. Mercator Map projection

EDIT Created a working example in PHP (because I suck at Java)

https://github.com/mfeldheim/mapStuff.git

EDIT2

Nice animation of the Mercator projection https://amp-reddit-com.cdn.ampproject.org/v/s/amp.reddit.com/r/educationalgifs/comments/5lhk8y/how_the_mercator_projection_distorts_the_poles/?usqp=mq331AQJCAEoAVgBgAEB&amp_js_v=0.1

like image 66
Michel Feldheim Avatar answered Sep 20 '22 07:09

Michel Feldheim


You cannot merely transpose from longitude/latitude to x/y like that because the world isn't flat. Have you look at this post? Converting longitude/latitude to X/Y coordinate

UPDATE - 1/18/13

I decided to give this a stab, and here's how I do it:-

public class MapService {     // CHANGE THIS: the output path of the image to be created     private static final String IMAGE_FILE_PATH = "/some/user/path/map.png";      // CHANGE THIS: image width in pixel     private static final int IMAGE_WIDTH_IN_PX = 300;      // CHANGE THIS: image height in pixel     private static final int IMAGE_HEIGHT_IN_PX = 500;      // CHANGE THIS: minimum padding in pixel     private static final int MINIMUM_IMAGE_PADDING_IN_PX = 50;      // formula for quarter PI     private final static double QUARTERPI = Math.PI / 4.0;      // some service that provides the county boundaries data in longitude and latitude     private CountyService countyService;      public void run() throws Exception {         // configuring the buffered image and graphics to draw the map         BufferedImage bufferedImage = new BufferedImage(IMAGE_WIDTH_IN_PX,                                                         IMAGE_HEIGHT_IN_PX,                                                         BufferedImage.TYPE_INT_RGB);          Graphics2D g = bufferedImage.createGraphics();         Map<RenderingHints.Key, Object> map = new HashMap<RenderingHints.Key, Object>();         map.put(RenderingHints.KEY_INTERPOLATION, RenderingHints.VALUE_INTERPOLATION_BICUBIC);         map.put(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);         map.put(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);         RenderingHints renderHints = new RenderingHints(map);         g.setRenderingHints(renderHints);          // min and max coordinates, used in the computation below         Point2D.Double minXY = new Point2D.Double(-1, -1);         Point2D.Double maxXY = new Point2D.Double(-1, -1);          // a list of counties where each county contains a list of coordinates that form the county boundary         Collection<Collection<Point2D.Double>> countyBoundaries = new ArrayList<Collection<Point2D.Double>>();          // for every county, convert the longitude/latitude to X/Y using Mercator projection formula         for (County county : countyService.getAllCounties()) {             Collection<Point2D.Double> lonLat = new ArrayList<Point2D.Double>();              for (CountyBoundary countyBoundary : county.getCountyBoundaries()) {                 // convert to radian                 double longitude = countyBoundary.getLongitude() * Math.PI / 180;                 double latitude = countyBoundary.getLatitude() * Math.PI / 180;                  Point2D.Double xy = new Point2D.Double();                 xy.x = longitude;                 xy.y = Math.log(Math.tan(QUARTERPI + 0.5 * latitude));                  // The reason we need to determine the min X and Y values is because in order to draw the map,                 // we need to offset the position so that there will be no negative X and Y values                 minXY.x = (minXY.x == -1) ? xy.x : Math.min(minXY.x, xy.x);                 minXY.y = (minXY.y == -1) ? xy.y : Math.min(minXY.y, xy.y);                  lonLat.add(xy);             }              countyBoundaries.add(lonLat);         }          // readjust coordinate to ensure there are no negative values         for (Collection<Point2D.Double> points : countyBoundaries) {             for (Point2D.Double point : points) {                 point.x = point.x - minXY.x;                 point.y = point.y - minXY.y;                  // now, we need to keep track the max X and Y values                 maxXY.x = (maxXY.x == -1) ? point.x : Math.max(maxXY.x, point.x);                 maxXY.y = (maxXY.y == -1) ? point.y : Math.max(maxXY.y, point.y);             }         }          int paddingBothSides = MINIMUM_IMAGE_PADDING_IN_PX * 2;          // the actual drawing space for the map on the image         int mapWidth = IMAGE_WIDTH_IN_PX - paddingBothSides;         int mapHeight = IMAGE_HEIGHT_IN_PX - paddingBothSides;          // determine the width and height ratio because we need to magnify the map to fit into the given image dimension         double mapWidthRatio = mapWidth / maxXY.x;         double mapHeightRatio = mapHeight / maxXY.y;          // using different ratios for width and height will cause the map to be stretched. So, we have to determine         // the global ratio that will perfectly fit into the given image dimension         double globalRatio = Math.min(mapWidthRatio, mapHeightRatio);          // now we need to readjust the padding to ensure the map is always drawn on the center of the given image dimension         double heightPadding = (IMAGE_HEIGHT_IN_PX - (globalRatio * maxXY.y)) / 2;         double widthPadding = (IMAGE_WIDTH_IN_PX - (globalRatio * maxXY.x)) / 2;          // for each country, draw the boundary using polygon         for (Collection<Point2D.Double> points : countyBoundaries) {             Polygon polygon = new Polygon();              for (Point2D.Double point : points) {                 int adjustedX = (int) (widthPadding + (point.getX() * globalRatio));                  // need to invert the Y since 0,0 starts at top left                 int adjustedY = (int) (IMAGE_HEIGHT_IN_PX - heightPadding - (point.getY() * globalRatio));                  polygon.addPoint(adjustedX, adjustedY);             }              g.drawPolygon(polygon);         }          // create the image file         ImageIO.write(bufferedImage, "PNG", new File(IMAGE_FILE_PATH));     } } 

RESULT: Image width = 600px, Image height = 600px, Image padding = 50px

enter image description here

RESULT: Image width = 300px, Image height = 500px, Image padding = 50px

enter image description here

like image 33
limc Avatar answered Sep 23 '22 07:09

limc