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 (Swissgrid, CH1903 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:
[sourcecode language=”python” gutter=”false”]
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][/sourcecode]