Conversion swisstopo, CH1903 (LV95/LV03) and WGS83

TL;DR; Summary and code

Pass in either Latitude/Longitude to wgs84_to_ch1903 (which by default converts to CH1903+) or “Rechtswert” (x) and “Hochwert” (y) to ch1903_to_wgs84 (which detects if its CH1903+ based on the length/value of the passed digits. The most recent version is available as a github gist from me .

Python

import numpy as np

def deci2sexa(angle):

    angle = np.asarray(angle)
    # Extract DMS
    degrees = angle.astype(int)
    minutes = (angle-degrees*60).astype(int)
    seconds = (((angle-degrees)*60)-minutes)*60
    # Result sexagesimal seconds
    return seconds + minutes * 60.0 + degrees * 3600.0

def wgs84_to_ch1903(lat, lon, plus=True):

    lat, lon = deci2sexa(lat), deci2sexa(lon)
    # Auxiliary values (% Bern)
    lat_aux = (lat - 169028.66) / 10000
    lng_aux = (lon - 26782.5) / 10000

    x = (200147.07 +
         308807.95 * lat_aux  +
         3745.25 * np.power(lng_aux, 2) +
         76.63 * np.power(lat_aux, 2) -
         194.56 * np.power(lng_aux, 2) * lat_aux +
         119.79 * np.power(lat_aux, 3))

    y = (600072.37 +
         211455.93 * lng_aux -
         10938.51 * lng_aux * lat_aux -
         0.36 * lng_aux * np.power(lat_aux, 2) -
         44.54 * np.power(lng_aux, 3))

    if plus:
        x += 1000000
        y += 2000000

    return x, y
    
  def ch1903_to_wgs84(x, y, plus='auto'):

    if plus == 'auto':
        if np.nanmax(x) > 1200000 or np.nanmax(y) > 2600000:
            plus = True
        else:
            plus = False

    #  Auxiliary values (% Bern)
    y_aux = (y - 600000)/1000000 # would be 2200000 for ch1903plus
    x_aux = (x - 200000)/1000000 # would be 1200000 for ch1903plus

    if plus:
        x_aux -= 1 # new ch1903plus system has another digit to distinguish it
        y_aux -= 2 # new ch1903plus system has another digit to distinguish it
    # Process lat
    lat = (16.9023892 +
           3.238272 * x_aux -
           0.270978 * np.power(y_aux, 2) -
           0.002528 * np.power(x_aux, 2) -
           0.0447 * np.power(y_aux, 2) * x_aux -
           0.0140 * np.power(x_aux, 3))

    # Process lng
    lon = (2.6779094 +
           4.728982 * y_aux +
           0.791484 * y_aux * x_aux +
           0.1306 * y_aux * np.power(x_aux, 2) -
           0.0436 * np.power(y_aux, 3))

	# Unit 10000" to 1 " and converts seconds to degrees (dec)
    lon = lon * 100 / 36
    lat = lat * 100 / 36
    return lat, lon

Quite often I find myself having to convert WGS84 to CH1903 coordinate systems and vice versa. Sometimes I even got neither and simply have a center point and some distance (looking at you ground-based remote sensing data). While swisstopo used to have (in 2017 or so) a library to download, the current easiest way is to actually use the github repo from Valentin Minder which contains converter for several programming languages.

However, the repo contains classes (which are great of course) but I often prefer a direct function (which in Python is also a class, but well …). As such, I used the same formulas you can find elsewhere from swisstopo to do the calculation. Since I usually do not care that much about the altitude in these cases there is no option (as of yet) to include it. One upside of this version of the conversion is that its array enabled, i.e. the respective coordinates can be either a single scalar or a numpy array (not a list though ;-)).

Conversion swisstopo, CH1903 (LV95/LV03) and WGS83

TL;DR; Summary and code

Further reading

One Reply to “Conversion swisstopo, CH1903 (LV95/LV03) and WGS83”