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

Question

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.

Cloots

Your formulas are just linear interpolation, which effectively implies you're doing an equirectangular projection rather than a Mercator. – Sesquiplane

I've updated the code, although still having problems with latitude – Cloots

As @Drew mentioned, if your map is a Marcator projection, you'll need to convert the lat/lng into x/y using a Mercator projection. Check if your map is Transverse Mercator or Spherical Mercator, then we'll get to the formulars... – Favian

It's a spherical Mercator projection – Cloots

Answer 1

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.]

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

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

Amor

SWEET iv been looking for this equation put simply in pseudo code for aaaages ty – Vitelline

@Michel: Can you add some example, please? – Cheerful

If i want "zoom" in my coordinate points, Where i can multiplicate the zoom in this algorithm? – Sebi

I don't see where that pseudo code takes into account the lat/long coordinates of the corners of the 2D map. Isn't that necessary? – Hominid

@MichelFeldheim Thanks so much for the github link. Future visitors, be warned you need the php-image-text package with php sudo apt-get install php-image-text in order to run this example. – Glutton

Is it on purpose that y depends also on the mapWidth value - in your pseudo code it says: y = (mapHeight/2)-(mapWidthmercN/(2*PI)); - shouldn't that be: y = (mapHeight/2)-(mapHeightmercN/(2*PI)); ? – Adulterer

Yes, the formula is simplified and not 100% accurate but it should be fine for common usage where you don't need scientific accuracy – Amor

What do we need to do make this work with a cropped map - the one whose bounds are not a full world map? – Guacin

Do you know how to convert the height value to appropriate z coordinate? – Backset

@Wojciech Danilo - you mean reverse calculation - pixel position to geographic lat/lon? – Amor

@MichelFeldheim: No, I ment something else. If ve've got (lat, lon, elevation), than we can convert it to (x,y,z), where z vould be matched against the local metrics of x and y. What I mean is that distance between two (x,y) (x+1,y+1) points after the projection in euclidean metrics depends strongly on their position. The same way the z value would vary, am I right? – Backset

Well - elevation is usually in metres above NN, so, depending on your visualisation it's a simple rule of three +/- offset using earth radius and pixel height of your landscape profile or 3d model – Amor

Longitude is spot on for me but the latitude is completely off. Mind you I'm trying to display world cities on a 158 x 68 array, maybe the small height explains the latitude offset. – Tiffanietiffanle

@Tiffanietiffanle please try this modified formula y = (mapHeight/2)-(mapHeight*mercN/(2*PI));. I think @Adulterer is actually right. The formula did work in my example, because I have used a map with same width and height – Amor

It's better! London and Paris appear exactly where I would expect them. Mexico City appears around Texas in the US and Buenos Aires in Venuzuela/Guyana. (Only tried those four cities so far). – Tiffanietiffanle

Success! I used your formula for the longitude, and for the latitude I used this simpler: y = h/2 -(latitude * h) / 180. It's approximate, but accurate enough for this application. – Tiffanietiffanle

and if you used this method to get an xy pair, how could you convert them back to LatLng? – Blasius

convert back is (javascript language): long = ((360 * x) / mapWidth) - 180; lat = 90 * (-1 + (4 * Math.atan(Math.pow(Math.E, (Math.PI - (2 * Math.PI * y) / mapHeight)))) / Math.PI); – Lupulin

What is the variable "R"? – Ieper

Its the radius of earth 6,371km. – Redoubtable

What is ln() in javascript? – Co

@Co Math.log() – Amor

Answer 2

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

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

Cardew

I've had a look, it doesn't seem like the author's code takes into consideration the size of the map image? I've updated my question with hopefully more accurate code-although still not right. – Cloots

I updated my post... this code takes account of the specified image width and height. – Cardew

Answer 3

Java version of original Google Maps JavaScript API v3 java script code is as following, it works with no problem

public final class GoogleMapsProjection2 
{
    private final int TILE_SIZE = 256;
    private PointF _pixelOrigin;
    private double _pixelsPerLonDegree;
    private double _pixelsPerLonRadian;

    public GoogleMapsProjection2()
    {
        this._pixelOrigin = new PointF(TILE_SIZE / 2.0,TILE_SIZE / 2.0);
        this._pixelsPerLonDegree = TILE_SIZE / 360.0;
        this._pixelsPerLonRadian = TILE_SIZE / (2 * Math.PI);
    }

    double bound(double val, double valMin, double valMax)
    {
        double res;
        res = Math.max(val, valMin);
        res = Math.min(res, valMax);
        return res;
    }

    double degreesToRadians(double deg) 
    {
        return deg * (Math.PI / 180);
    }

    double radiansToDegrees(double rad) 
    {
        return rad / (Math.PI / 180);
    }

    PointF fromLatLngToPoint(double lat, double lng, int zoom)
    {
        PointF point = new PointF(0, 0);

        point.x = _pixelOrigin.x + lng * _pixelsPerLonDegree;       

        // Truncating to 0.9999 effectively limits latitude to 89.189. This is
        // about a third of a tile past the edge of the world tile.
        double siny = bound(Math.sin(degreesToRadians(lat)), -0.9999,0.9999);
        point.y = _pixelOrigin.y + 0.5 * Math.log((1 + siny) / (1 - siny)) *- _pixelsPerLonRadian;

        int numTiles = 1 << zoom;
        point.x = point.x * numTiles;
        point.y = point.y * numTiles;
        return point;
     }

    PointF fromPointToLatLng(PointF point, int zoom)
    {
        int numTiles = 1 << zoom;
        point.x = point.x / numTiles;
        point.y = point.y / numTiles;       

        double lng = (point.x - _pixelOrigin.x) / _pixelsPerLonDegree;
        double latRadians = (point.y - _pixelOrigin.y) / - _pixelsPerLonRadian;
        double lat = radiansToDegrees(2 * Math.atan(Math.exp(latRadians)) - Math.PI / 2);
        return new PointF(lat, lng);
    }

    public static void main(String []args) 
    {
        GoogleMapsProjection2 gmap2 = new GoogleMapsProjection2();

        PointF point1 = gmap2.fromLatLngToPoint(41.850033, -87.6500523, 15);
        System.out.println(point1.x+"   "+point1.y);
        PointF point2 = gmap2.fromPointToLatLng(point1,15);
        System.out.println(point2.x+"   "+point2.y);
    }
}

public final class PointF 
{
    public double x;
    public double y;

    public PointF(double x, double y)
    {
        this.x = x;
        this.y = y;
    }
}

Erskine

What defines the tile size? I am trying to use something like this in Unity. Is it the image size? or is it a constant for google maps or mercator projection on the web? – Incur

tile size is the image size of map tile, when you receive a map tile image from google servers it is a 256x256 size of png file – Erskine

This works for latitude, but longitude values is way off. What is missing? – Bramblett

Answer 4

JAVA only?

Python code here! Refer to Convert latitude/longitude point to a pixels (x,y) on mercator projection

import math
from numpy import log as ln

# Define the size of map
mapWidth    = 200
mapHeight   = 100


def convert(latitude, longitude):
    # get x value
    x = (longitude + 180) * (mapWidth / 360)

    # convert from degrees to radians
    latRad = (latitude * math.pi) / 180

    # get y value
    mercN = ln(math.tan((math.pi / 4) + (latRad / 2)))
    y     = (mapHeight / 2) - (mapWidth * mercN / (2 * math.pi))
    
    return x, y

print(convert(41.145556, 121.2322))

Answer:

(167.35122222222225, 24.877939817552335)

Boughten

Answer 5

 public static String getTileNumber(final double lat, final double lon, final int zoom) {
 int xtile = (int)Math.floor( (lon + 180) / 360 * (1<<zoom) ) ;
 int ytile = (int)Math.floor( (1 - Math.log(Math.tan(Math.toRadians(lat)) + 1 /  Math.cos(Math.toRadians(lat))) / Math.PI) / 2 * (1<<zoom) ) ;
if (xtile < 0)
 xtile=0;
if (xtile >= (1<<zoom))
 xtile=((1<<zoom)-1);
if (ytile < 0)
 ytile=0;
if (ytile >= (1<<zoom))
 ytile=((1<<zoom)-1);
return("" + zoom + "/" + xtile + "/" + ytile);
 }
}

Soil

Answer 6

I'm new here, just to write, as I've been following the community for some years. I'm happy to be able to contribute.

Well, it took me practically a day in search of that and your question encouraged me to continue the search.

I arrived at the following function, which works! Credits for this article: https://towardsdatascience.com/geotiff-coordinate-querying-with-javascript-5e6caaaf88cf

var bbox = [minLong, minLat, maxLong, maxLat];
var pixelWidth = mapWidth;
var pixelHeight = mapHeight;
var bboxWidth = bbox[2] - bbox[0];
var bboxHeight = bbox[3] - bbox[1];

var convertToXY = function(latitude, longitude) {
    var widthPct = ( longitude - bbox[0] ) / bboxWidth;
    var heightPct = ( latitude - bbox[1] ) / bboxHeight;
    var x = Math.floor( pixelWidth * widthPct );
    var y = Math.floor( pixelHeight * ( 1 - heightPct ) );
    return { x, y };
}

Polyphagia