2

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.
    deg_to_rad = math.pi / 180.0
    phi = latitude * deg_to_rad
    lambda_ = longitude * deg_to_rad
    lambda_zero = central_meridian * deg_to_rad

    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.
    deg_to_rad = math.pi / 180
    lambda_zero = central_meridian * deg_to_rad
    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))
    lon_radian = lambda_zero + delta_lambda
    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.

  • 3
    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 – Spacedman Aug 13 '18 at 18:22
  • That JS seems to have options for lots of SWEREF99 projections - are they not what you have? gist.github.com/plopp/… – Spacedman Aug 13 '18 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 Aug 13 '18 at 18:45
  • 2
    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. – mkennedy Aug 13 '18 at 20:07
  • 1
    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. – mkennedy Aug 14 '18 at 0:48
2

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.
    deg_to_rad = math.pi / 180.0
    phi = latitude * deg_to_rad
    lambda = longitude * deg_to_rad

    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 Aug 14 '18 at 17:30

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.