Projection from WGS1984 to Swiss national grid (CH1903 LV03)

In my last project I used (standalone) Python for geoprocessing. Since ArcGIS or something like that was not there to help with projecting geodata from one coordinate system to another, I wrote a function which converts well-known WGS1984 to Swiss national grid coordinates (SwissgridCH1903 LV03).

Swisstopo has the formulas of an approximate transformation and also funtions in C#, Javascript and PHP, but not in Python. So here is to your disposal:

def projWGS1984ToCH1903(x, y):
  """Converts coordinates from WGS1984 to CH1903_LV03 using the
     method from http://www.swisstopo.admin.ch/internet/swisstopo/
     de/home/products/software/products/skripts.html, returns a
     list [x,y]
     :param x: x coordinates in degrees in WGS1984 :
     :param y: y coordinates in degrees in WGS1984 """

  # Transformation into sexagesimal seconds
  x = x * 3600
  y = y * 3600

  # Latitude and longitude difference to Bern
  x_fact = (x - 26782.5) / 10000 # LAMBDA
  y_fact = (y - 169028.66) / 10000 # PHI

  x = 600072.37 + 211455.93 * x_fact - 10938.51 * x_fact * y_fact - 0.36 * x_fact * y_fact**2 - 44.54 * x_fact**3
  y = 200147.07 + 308807.95 * y_fact + 3745.25 * x_fact**2 + 76.63 * y_fact**2 - 194.56 * x_fact**2 * y_fact + 119.79 * y_fact**3

  return [x, y]
Advertisements

Leave a reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s