Geographic Information Systems Asked by user50614 on July 17, 2021
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?
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;
}
Answered by Russell at ISC on July 17, 2021
from pyproj import Proj, transform
P3857 = Proj(init='epsg:3857')
P4326 = Proj(init='epsg:4326')
x,y = transform(P4326, P3857, lon, lat)
Answered by Marcel Wilson on July 17, 2021
Following up on the answer of @Russel at ICS. 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
}
Answered by Wilmar van Ommeren on July 17, 2021
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
Answered by Hua Zhang on July 17, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP