Converting latitude, longitude (EPSG:4326) into EPSG:3857?

Question

I'm looking to study the math involved to see how a given latitude and longitude in decimal degrees can be converted to EPSG:3857. Can anybody point to a good reference or a possible Ppython/C/C++/Java library source code for checking the math?

Answer 1

from pyproj import Proj, transform

P3857 = Proj(init='epsg:3857')
P4326 = Proj(init='epsg:4326')

x,y = transform(P4326, P3857, lon, lat)

Answer 2

This code snippet is C# and using an array called vertex to hold [x, y]

double smRadius = 6378136.98;
double smRange = smRadius * Math.PI * 2.0;
double smLonToX = smRange / 360.0;
double smRadiansOverDegrees = Math.PI / 180.0;

...

// compute x-map-unit
vertex[0] *= smLonToX;

double y = vertex[1];

// compute y-map-unit
if (y > 86.0)
{
    vertex[1] = smRange;
}
else if (y < -86.0)
{
    vertex[1] = -smRange;
}
else
{
    y *= smRadiansOverDegrees;
    y = Math.Log(Math.Tan(y) + (1.0 / Math.Cos(y)), Math.E);
    vertex[1] = y * smRadius; 
}

Answer 3

Here is my Python code. Ref to this link

#   Converting lat, lon (epsg:4326) into EPSG:3857
# Ref: https://stackoverflow.com/questions/37523872/converting-coordinates-from-epsg-3857-to-4326-dotspatial/40403522#40403522

def ConvertCoordinate(lon, lat):
    lonInEPSG4326 = lon
    latInEPSG4326 = lat

    lonInEPSG3857 = (lonInEPSG4326 * 20037508.34 / 180)
    latInEPSG3857 = (math.log(math.tan((90 + latInEPSG4326) * math.pi / 360)) / (math.pi / 180)) * (20037508.34 / 180)

    print("{0},{1}".format(lonInEPSG3857,latInEPSG3857))
    return lonInEPSG3857, latInEPSG3857

Answer 4

I use the osgro library, using the ogr and osr classes. This way you can transform anything to anything (pretty much)

from osgeo import ogr, osr

def transform_4326_to_3857(lat, long):
    # make spatial references
    # 4326 in
    InSR = osr.SpatialReference()
    InSR.ImportFromEPSG(4326)
    # 3857 out
    OutSR = osr.SpatialReference()
    OutSR.ImportFromEPSG(3857)
    
    # create the point
    Point = ogr.Geometry(ogr.wkbPoint)
    Point.AddPoint(lat,long)

    # assign the spacial reference to the point
    Point.AssignSpatialReference(InSR)

    # transform
    transform_error_code = Point.TransformTo(OutSR)

    # error check (eg, point is out of bounds, etc)
    if transform_error_code == 0:
        # success
        return (Point.GetX(),Point.GetY())
    else:
        #failed
        logging.warn("Transform failed! Error {transform_error_code}")
        return None

Answer 5

Following up on the answer of @Russell at ISC. Today I rewrote his code to convert WGS84 (EPSG:4326) coordinates to WGS84/Pseudo Mercator (EPSG:3857) in R:

vertex = list(x_coordinate, y_coordinate)
smRadius = 6378136.98
smRange = smRadius * pi * 2.0
smLonToX = smRange / 360.0
smRadiansOverDegrees = pi / 180.0

vertex[[1]] = vertex[[1]] *smLonToX

y = vertex[[2]]

if (y > 86.0){
  vertex[[2]] = smRange
} else if (y < -86.0){
  vertex[[2]] = -smRange
} else {
  y = y * smRadiansOverDegrees
  y = log(tan(y) + (1.0 / cos(y)), base = exp(1))
  vertex[[2]] = y * smRadius
}