12

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?

2

5 Answers 5

6
from pyproj import Proj, transform

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

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

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

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
0

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
0

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
}

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.