From SWEREF99 to EPSG:3857

I have a huge dataset of coordinates referenced with the SWEREF99 coordinate system. Now I would like to convert all these coordinates (and other data as well) to a coordinate system that can be rendered easily with OpenStreetMap/Leaflet. For that I have a python script that goes through the csvs, but I just can't figure out the formula that would help me convert to EPSG:3857.

I have found a code (in Javascript) that converts to RT 90 from Arnold Andreasson. I have adapted this code to python (somehow) but even from there after hours of browsing through here and on google I can't seem to find something that would convert from RT90 or SWEREF99 to EPSG:3857.

Small information: using an online service is definitely not an option, sensitive data + huge dataset scattered in several CSVs here.

I have seen that some answers here rely on proj4 but I've had trouble finding a good python version of that that can be installed on my Mac... Also, I'd like to rely on external libraries as less as possible.

Here is the python code:

``````import math

axis = None # Semi-major axis of the ellipsoid.
flattening = None # Flattening of the ellipsoid.
central_meridian = None # Central meridian for the projection.
lat_of_origin = None # Latitude of origin.
scale = None # Scale on central meridian.
false_northing = None # Offset for origo.
false_easting = None # Offset for origo.

# Parameters for RT90 and SWEREF99TM.
# Note: Parameters for RT90 are choosen to eliminate the
# differences between Bessel and GRS80-ellipsoides.
# Bessel-iants should only be used if lat/long are given as
# RT90-lat/long based on the Bessel ellipsoide (from old maps).
# Parameter: projection (string). Must match if-statement.
def swedish_params(projection) :

global central_meridian
global scale
global false_northing
global false_easting
global lat_of_origin

# RT90 parameters, GRS 80 ellipsoid.
if (projection == "rt90_7.5_gon_v") :
grs80_params()
central_meridian = 11.0 + 18.375/60.0
scale = 1.000006000000
false_northing = -667.282
false_easting = 1500025.141

elif (projection == "rt90_5.0_gon_v") :
grs80_params()
central_meridian = 13.0 + 33.376/60.0
scale = 1.000005800000
false_northing = -667.130
false_easting = 1500044.695

elif (projection == "rt90_2.5_gon_v") :
grs80_params()
central_meridian = 15.0 + 48.0/60.0 + 22.624306/3600.0
scale = 1.00000561024
false_northing = -667.711
false_easting = 1500064.274

elif (projection == "rt90_0.0_gon_v") :
grs80_params()
central_meridian = 18.0 + 3.378/60.0
scale = 1.000005400000
false_northing = -668.844
false_easting = 1500083.521

elif (projection == "rt90_2.5_gon_o") :
grs80_params()
central_meridian = 20.0 + 18.379/60.0
scale = 1.000005200000
false_northing = -670.706
false_easting = 1500102.765

elif (projection == "rt90_5.0_gon_o") :
grs80_params()
central_meridian = 22.0 + 33.380/60.0
scale = 1.000004900000
false_northing = -672.557
false_easting = 1500121.846

# RT90 parameters, Bessel 1841 ellipsoid.
elif (projection == "bessel_rt90_7.5_gon_v") :
bessel_params()
central_meridian = 11.0 + 18.0/60.0 + 29.8/3600.0

elif (projection == "bessel_rt90_5.0_gon_v") :
bessel_params()
central_meridian = 13.0 + 33.0/60.0 + 29.8/3600.0

elif (projection == "bessel_rt90_2.5_gon_v") :
bessel_params()
central_meridian = 15.0 + 48.0/60.0 + 29.8/3600.0

elif (projection == "bessel_rt90_0.0_gon_v") :
bessel_params()
central_meridian = 18.0 + 3.0/60.0 + 29.8/3600.0

elif (projection == "bessel_rt90_2.5_gon_o") :
bessel_params()
central_meridian = 20.0 + 18.0/60.0 + 29.8/3600.0

elif (projection == "bessel_rt90_5.0_gon_o") :
bessel_params()
central_meridian = 22.0 + 33.0/60.0 + 29.8/3600.0

# SWEREF99TM and SWEREF99ddmm  parameters.
elif (projection == "sweref_99_tm") :
sweref99_params()
central_meridian = 15.00
lat_of_origin = 0.0
scale = 0.9996
false_northing = 0.0
false_easting = 500000.0

elif (projection == "sweref_99_1200") :
sweref99_params()
central_meridian = 12.00

elif (projection == "sweref_99_1330") :
sweref99_params()
central_meridian = 13.50

elif (projection == "sweref_99_1500") :
sweref99_params()
central_meridian = 15.00

elif (projection == "sweref_99_1630") :
sweref99_params()
central_meridian = 16.50

elif (projection == "sweref_99_1800") :
sweref99_params()
central_meridian = 18.00

elif (projection == "sweref_99_1415") :
sweref99_params()
central_meridian = 14.25

elif (projection == "sweref_99_1545") :
sweref99_params()
central_meridian = 15.75

elif (projection == "sweref_99_1715") :
sweref99_params()
central_meridian = 17.25

elif (projection == "sweref_99_1845") :
sweref99_params()
central_meridian = 18.75

elif (projection == "sweref_99_2015") :
sweref99_params()
central_meridian = 20.25

elif (projection == "sweref_99_2145") :
sweref99_params()
central_meridian = 21.75

elif (projection == "sweref_99_2315") :
sweref99_params()
central_meridian = 23.25

# Test-case:
#  Lat: 66 0'0", lon: 24 0'0".
#  X:1135809.413803 Y:555304.016555.
elif (projection == "test_case") :
axis = 6378137.0
flattening = 1.0 / 298.257222101
central_meridian = 13.0 + 35.0/60.0 + 7.692000/3600.0
lat_of_origin = 0.0
scale = 1.000002540000
false_northing = -6226307.8640
false_easting = 84182.8790

# Not a valid projection.
else :
central_meridian = None

# Sets of default parameters.
def grs80_params() :

global axis
global flattening
global central_meridian
global lat_of_origin

axis = 6378137.0 # GRS 80.
flattening = 1.0 / 298.257222101 # GRS 80.
central_meridian = None
lat_of_origin = 0.0

def bessel_params() :

global axis
global flattening
global central_meridian
global lat_of_origin
global scale
global false_northing
global false_easting

axis = 6377397.155 # Bessel 1841.
flattening = 1.0 / 299.1528128 # Bessel 1841.
central_meridian = None
lat_of_origin = 0.0
scale = 1.0
false_northing = 0.0
false_easting = 1500000.0

def sweref99_params() :

global axis
global flattening
global central_meridian
global lat_of_origin
global scale
global false_northing
global false_easting

axis = 6378137.0 # GRS 80.
flattening = 1.0 / 298.257222101 # GRS 80.
central_meridian = None
lat_of_origin = 0.0
scale = 1.0
false_northing = 0.0
false_easting = 150000.0

# Conversion from geodetic coordinates to grid coordinates.
def geodetic_to_grid(latitude, longitude) :
x_y = [0] * 2
if (central_meridian == None) :
return x_y

# Prepare ellipsoid-based stuff.
e2 = flattening * (2.0 - flattening)
n = flattening / (2.0 - flattening)
a_roof = axis / (1.0 + n) * (1.0 + n*n/4.0 + n*n*n*n/64.0)
A = e2
B = (5.0*e2*e2 - e2*e2*e2) / 6.0
C = (104.0*e2*e2*e2 - 45.0*e2*e2*e2*e2) / 120.0
D = (1237.0*e2*e2*e2*e2) / 1260.0
beta1 = n/2.0 - 2.0*n*n/3.0 + 5.0*n*n*n/16.0 + 41.0*n*n*n*n/180.0
beta2 = 13.0*n*n/48.0 - 3.0*n*n*n/5.0 + 557.0*n*n*n*n/1440.0
beta3 = 61.0*n*n*n/240.0 - 103.0*n*n*n*n/140.0
beta4 = 49561.0*n*n*n*n/161280.0

# Convert.

phi_star = phi - math.sin(phi) * math.cos(phi) * (A + B*math.pow(math.sin(phi), 2) + C*math.pow(math.sin(phi), 4) + D*math.pow(math.sin(phi), 6))
delta_lambda = lambda_ - lambda_zero
xi_prim = math.atan(math.tan(phi_star) / math.cos(delta_lambda))
eta_prim = math_atanh(math.cos(phi_star) * math.sin(delta_lambda))
x = scale * a_roof * (xi_prim +beta1 * math.sin(2.0*xi_prim) * math_cosh(2.0*eta_prim) +beta2 * math.sin(4.0*xi_prim) * math_cosh(4.0*eta_prim) +beta3 * math.sin(6.0*xi_prim) * math_cosh(6.0*eta_prim) +beta4 * math.sin(8.0*xi_prim) * math_cosh(8.0*eta_prim)) + false_northing
y = scale * a_roof * (eta_prim +beta1 * math.cos(2.0*xi_prim) * math_sinh(2.0*eta_prim) +beta2 * math.cos(4.0*xi_prim) * math_sinh(4.0*eta_prim) +beta3 * math.cos(6.0*xi_prim) * math_sinh(6.0*eta_prim) +beta4 * math.cos(8.0*xi_prim) * math_sinh(8.0*eta_prim)) + false_easting
x_y[0] = math.round(x * 1000.0) / 1000.0
x_y[1] = math.round(y * 1000.0) / 1000.0
#  x_y[0] = x
#  x_y[1] = y
return x_y

# Conversion from grid coordinates to geodetic coordinates.
def grid_to_geodetic(x, y) :
lat_lon = [0] * 2
if (central_meridian == None) :
return lat_lon

# Prepare ellipsoid-based stuff.
e2 = flattening * (2.0 - flattening)
n = flattening / (2.0 - flattening)
a_roof = axis / (1.0 + n) * (1.0 + n*n/4.0 + n*n*n*n/64.0)
delta1 = n/2.0 - 2.0*n*n/3.0 + 37.0*n*n*n/96.0 - n*n*n*n/360.0
delta2 = n*n/48.0 + n*n*n/15.0 - 437.0*n*n*n*n/1440.0
delta3 = 17.0*n*n*n/480.0 - 37*n*n*n*n/840.0
delta4 = 4397.0*n*n*n*n/161280.0

Astar = e2 + e2*e2 + e2*e2*e2 + e2*e2*e2*e2
Bstar = -(7.0*e2*e2 + 17.0*e2*e2*e2 + 30.0*e2*e2*e2*e2) / 6.0
Cstar = (224.0*e2*e2*e2 + 889.0*e2*e2*e2*e2) / 120.0
Dstar = -(4279.0*e2*e2*e2*e2) / 1260.0

# Convert.
xi = (x - false_northing) / (scale * a_roof)
eta = (y - false_easting) / (scale * a_roof)
xi_prim = xi - delta1*math.sin(2.0*xi) * math_cosh(2.0*eta) - delta2*math.sin(4.0*xi) * math_cosh(4.0*eta) - delta3*math.sin(6.0*xi) * math_cosh(6.0*eta) - delta4*math.sin(8.0*xi) * math_cosh(8.0*eta)
eta_prim = eta - delta1*math.cos(2.0*xi) * math_sinh(2.0*eta) - delta2*math.cos(4.0*xi) * math_sinh(4.0*eta) - delta3*math.cos(6.0*xi) * math_sinh(6.0*eta) - delta4*math.cos(8.0*xi) * math_sinh(8.0*eta)
phi_star = math.asin(math.sin(xi_prim) / math_cosh(eta_prim))
delta_lambda = math.atan(math_sinh(eta_prim) / math.cos(xi_prim))
lat_radian = phi_star + math.sin(phi_star) * math.cos(phi_star) * (Astar + Bstar*math.pow(math.sin(phi_star), 2) + Cstar*math.pow(math.sin(phi_star), 4) + Dstar*math.pow(math.sin(phi_star), 6))
lat_lon[0] = lat_radian * 180.0 / math.pi
lat_lon[1] = lon_radian * 180.0 / math.pi
return lat_lon

# Missing defs in the math library.
def math_sinh(value) :
return 0.5 * (math.exp(value) - math.exp(-value))

def math_cosh(value) :
return 0.5 * (math.exp(value) + math.exp(-value))

def math_atanh(value) :
return 0.5 * math.log((1.0 + value) / (1.0 - value))
``````

I unfortunately do not have more information about the coordinate system... All I know is what I told you: SWEREF99. But I suppose it must be the most standard one (is there?). I'm new to all these coordinate systems so I do apologise if I do not provide all the necessary information at ounce.

• Plenty of people have no problem with Proj.4 on their Macs. Why not fix that problem? If you could use Proj4, is this the specification of the projection you have: spatialreference.org/ref/epsg/3006 Commented Aug 13, 2018 at 18:22
• That JS seems to have options for lots of SWEREF99 projections - are they not what you have? gist.github.com/plopp/… Commented Aug 13, 2018 at 18:26
• @Spacedman Proj4.py doesn't work for me and as I stated I'd like not to rely on external libraries... Also, yes that code converts to RT90, not to EPSG:3857
– LBes
Commented Aug 13, 2018 at 18:45
• Is it SWEREF99 TM (AKA EPSG::3006)? Or one of the 20 or so other larger-scale projected CRS? Do code this yourself, you'll need to support transverse Mercator on an ellipsoid and Mercator on a sphere. SWEREF99 ~=WGS84 so you can skip the GeoCRS/datum transformation algorithms. Commented Aug 13, 2018 at 20:07
• All the RT90/SWEREF99 projected CRS use transverse Mercator so that must be what's implemented by Andreasson. EPSG:3857 uses Mercator on a sphere, see Snyder PDF page 41. Do you really need it? Google will take lat/lon values. Commented Aug 14, 2018 at 0:48

If the data is SWEREF99 TM, in Stockholm, the coordinates will be in the range of:

``````easting: 674000 m
northing: 6580000 m
``````

which is listed as "sweref_99_tm" in the code.

If it's in SWEREF99 18 00:

``````easting: 154000 m
northing: 6580000 m
``````

which is listed as "sweref_99_1800" in the code.

If it's in SWEREF99 0 gon emulation (of RT90 0 gon):

``````easting: 1500000 m
northing: 6580000 m
``````

which is listed as "rt90_0.0_gon_v" in the code.

Depending on which one you have, you would set the appropriate parameters and run the coordinates through grid_to_geodetic. You'll end up with latitude and longitude values in SWEREF99 (a geographic coordinate reference system). You can possibly stop here because latitude-longitude values should be able to be used in leaflet. I would call/reference them as WGS84.

If you still want to convert them to Web Mercator, EPSG::3857, you'll need to run them through geodetic_to_grid, but that code needs to change to be spherical Mercator.

``````# Conversion from geodetic coordinates to grid coordinates.
def geodetic_to_grid(latitude, longitude) :
x_y = [0] * 2
if (central_meridian == None) :
return x_y

# Prepare ellipsoid-based stuff.
axis = 6378137.0

# Convert.

x = axis * lambda

y = axis * math.log(math.tan(math.pi/4.0 + phi/2.0))

# I don't know if you need these
x_y[0] = math.round(x * 1000.0) / 1000.0
x_y[1] = math.round(y * 1000.0) / 1000.0

return x_y
``````

I can simplify the Mercator math quite a bit because 3857 has these parameter values:

``````central meridian: 0
false easting: 0
false northing: 0
latitude of true scale / standard parallel: 0
``````

Note: I don't know python so won't be any help if this doesn't work.

• Thank you very much for that answer. I will try this out and let you know. :)
– LBes
Commented Aug 14, 2018 at 17:30