185 lines
6.9 KiB
Python
185 lines
6.9 KiB
Python
#!/usr/bin/python2
|
|
#-*- coding: utf-8 -*-
|
|
|
|
# The MIT License (MIT)
|
|
#
|
|
# Copyright (c) 2014 Federal Office of Topography swisstopo, Wabern, CH and Aaron Schmocker
|
|
#
|
|
# Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
# of this software and associated documentation files (the "Software"), to deal
|
|
# in the Software without restriction, including without limitation the rights
|
|
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
# copies of the Software, and to permit persons to whom the Software is
|
|
# furnished to do so, subject to the following conditions:
|
|
#
|
|
# The above copyright notice and this permission notice shall be included in
|
|
# all copies or substantial portions of the Software.
|
|
#
|
|
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
|
# THE SOFTWARE.
|
|
#
|
|
|
|
# WGS84 <-> LV03 converter based on the scripts of swisstopo written for python2.7
|
|
# Aaron Schmocker [aaron@duckpond.ch]
|
|
# vim: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
|
|
|
|
# Source: http://www.mont-terri.ch/internet/swisstopo/en/home/products/software/products/skripts.html (see PDFs under "Documentation")
|
|
# Updated 10.10.2016
|
|
# Please validate your results with NAVREF on-line service: http://www.swisstopo.admin.ch/internet/swisstopo/en/home/apps/calc/navref.html (difference ~ 1-2m)
|
|
|
|
import math
|
|
|
|
class GPSConverter(object):
|
|
'''
|
|
GPS Converter class which is able to perform convertions between the
|
|
CH1903 and WGS84 system.
|
|
'''
|
|
# Convert CH y/x/h to WGS height
|
|
def CHtoWGSheight(self, y, x, h):
|
|
# Axiliary values (% Bern)
|
|
y_aux = (y - 600000) / 1000000
|
|
x_aux = (x - 200000) / 1000000
|
|
h = (h + 49.55) - (12.60 * y_aux) - (22.64 * x_aux)
|
|
return h
|
|
|
|
# Convert CH y/x to WGS lat
|
|
def CHtoWGSlat(self, y, x):
|
|
# Axiliary values (% Bern)
|
|
y_aux = (y - 600000) / 1000000
|
|
x_aux = (x - 200000) / 1000000
|
|
lat = (16.9023892 + (3.238272 * x_aux)) + \
|
|
- (0.270978 * pow(y_aux, 2)) + \
|
|
- (0.002528 * pow(x_aux, 2)) + \
|
|
- (0.0447 * pow(y_aux, 2) * x_aux) + \
|
|
- (0.0140 * pow(x_aux, 3))
|
|
# Unit 10000" to 1" and convert seconds to degrees (dec)
|
|
lat = (lat * 100) / 36
|
|
return lat
|
|
|
|
# Convert CH y/x to WGS long
|
|
def CHtoWGSlng(self, y, x):
|
|
# Axiliary values (% Bern)
|
|
y_aux = (y - 600000) / 1000000
|
|
x_aux = (x - 200000) / 1000000
|
|
lng = (2.6779094 + (4.728982 * y_aux) + \
|
|
+ (0.791484 * y_aux * x_aux) + \
|
|
+ (0.1306 * y_aux * pow(x_aux, 2))) + \
|
|
- (0.0436 * pow(y_aux, 3))
|
|
# Unit 10000" to 1" and convert seconds to degrees (dec)
|
|
lng = (lng * 100) / 36
|
|
return lng
|
|
|
|
# Convert decimal angle (° dec) to sexagesimal angle (dd.mmss,ss)
|
|
def DecToSexAngle(self, dec):
|
|
degree = int(math.floor(dec))
|
|
minute = int(math.floor((dec - degree) * 60))
|
|
second = (((dec - degree) * 60) - minute) * 60
|
|
return degree + (float(minute) / 100) + (second / 10000)
|
|
|
|
# Convert sexagesimal angle (dd.mmss,ss) to seconds
|
|
def SexAngleToSeconds(self, dms):
|
|
degree = 0
|
|
minute = 0
|
|
second = 0
|
|
degree = math.floor(dms)
|
|
minute = math.floor((dms - degree) * 100)
|
|
second = (((dms - degree) * 100) - minute) * 100
|
|
return second + (minute * 60) + (degree * 3600)
|
|
|
|
# Convert sexagesimal angle (dd.mmss) to decimal angle (degrees)
|
|
def SexToDecAngle(self, dms):
|
|
degree = 0
|
|
minute = 0
|
|
second = 0
|
|
degree = math.floor(dms)
|
|
minute = math.floor((dms - degree) * 100)
|
|
second = (((dms - degree) * 100) - minute) * 100
|
|
return degree + (minute / 60) + (second / 3600)
|
|
|
|
# Convert WGS lat/long (° dec) and height to CH h
|
|
def WGStoCHh(self, lat, lng, h):
|
|
lat = self.DecToSexAngle(lat)
|
|
lng = self.DecToSexAngle(lng)
|
|
lat = self.SexAngleToSeconds(lat)
|
|
lng = self.SexAngleToSeconds(lng)
|
|
# Axiliary values (% Bern)
|
|
lat_aux = (lat - 169028.66) / 10000
|
|
lng_aux = (lng - 26782.5) / 10000
|
|
h = (h - 49.55) + (2.73 * lng_aux) + (6.94 * lat_aux)
|
|
return h
|
|
|
|
# Convert WGS lat/long (° dec) to CH x
|
|
def WGStoCHx(self, lat, lng):
|
|
lat = self.DecToSexAngle(lat)
|
|
lng = self.DecToSexAngle(lng)
|
|
lat = self.SexAngleToSeconds(lat)
|
|
lng = self.SexAngleToSeconds(lng)
|
|
# Axiliary values (% Bern)
|
|
lat_aux = (lat - 169028.66) / 10000
|
|
lng_aux = (lng - 26782.5) / 10000
|
|
x = ((200147.07 + (308807.95 * lat_aux) + \
|
|
+ (3745.25 * pow(lng_aux, 2)) + \
|
|
+ (76.63 * pow(lat_aux,2))) + \
|
|
- (194.56 * pow(lng_aux, 2) * lat_aux)) + \
|
|
+ (119.79 * pow(lat_aux, 3))
|
|
return x
|
|
|
|
# Convert WGS lat/long (° dec) to CH y
|
|
def WGStoCHy(self, lat, lng):
|
|
lat = self.DecToSexAngle(lat)
|
|
lng = self.DecToSexAngle(lng)
|
|
lat = self.SexAngleToSeconds(lat)
|
|
lng = self.SexAngleToSeconds(lng)
|
|
# Axiliary values (% Bern)
|
|
lat_aux = (lat - 169028.66) / 10000
|
|
lng_aux = (lng - 26782.5) / 10000
|
|
y = (600072.37 + (211455.93 * lng_aux)) + \
|
|
- (10938.51 * lng_aux * lat_aux) + \
|
|
- (0.36 * lng_aux * pow(lat_aux, 2)) + \
|
|
- (44.54 * pow(lng_aux, 3))
|
|
return y
|
|
|
|
def LV03toWGS84(self, east, north, height):
|
|
'''
|
|
Convert LV03 to WGS84 Return a array of double that contain lat, long,
|
|
and height
|
|
'''
|
|
d = []
|
|
d.append(self.CHtoWGSlat(east, north))
|
|
d.append(self.CHtoWGSlng(east, north))
|
|
d.append(self.CHtoWGSheight(east, north, height))
|
|
return d
|
|
|
|
def WGS84toLV03(self, latitude, longitude, ellHeight):
|
|
'''
|
|
Convert WGS84 to LV03 Return an array of double that contaign east,
|
|
north, and height
|
|
'''
|
|
d = []
|
|
d.append(self.WGStoCHy(latitude, longitude))
|
|
d.append(self.WGStoCHx(latitude, longitude))
|
|
d.append(self.WGStoCHh(latitude, longitude, ellHeight))
|
|
return d
|
|
|
|
if __name__ == "__main__":
|
|
''' Example usage for the GPSConverter class.'''
|
|
|
|
converter = GPSConverter()
|
|
|
|
# Coordinates
|
|
wgs84 = [46.95108, 7.438637, 0]
|
|
lv03 = []
|
|
|
|
# Convert WGS84 to LV03 coordinates
|
|
lv03 = converter.WGS84toLV03(wgs84[0], wgs84[1], wgs84[2])
|
|
|
|
print "WGS84: "
|
|
print wgs84
|
|
print "LV03: "
|
|
print lv03
|