# -*- coding: utf-8 -*-
"""Utility methods to read GPS data and convert them via pyproj."""
import sys
import os
from math import floor
import numpy as np
import pygimli as pg
[docs]
def findUTMZone(lon, lat):
"""Find utm zone for lon and lat values.
lon -180 -- -174 -> 1 ... 174 -- 180 -> 60
lat < 0 hemisphere = S, > 0 hemisphere = N
Parameters
----------
lon : float
lat : float
Returns
-------
str :
zone + hemisphere
"""
zone = int((int(lon) + 180) / 6) + 1
if lat > 0:
return str(zone) + 'N'
return str(zone) + 'S'
[docs]
def getUTMProjection(zone, ellps='WGS84'):
"""Create and return the current coordinate projection.
This is a proxy for pyproj.
Parameters
----------
utmZone : str
Zone for for UTM
ellipsoid : str
ellipsoid based on ['wgs84']
Returns
-------
Pyproj Projection
"""
pyproj = pg.optImport('pyproj', 'Coordinate transformations.')
return pyproj.Proj(proj='utm', zone=zone, ellps=ellps)
[docs]
def getProjection(name, ref=None, **kwargs):
"""Syntactic sugar to get some default Projections."""
pyproj = pg.optImport('pyproj', 'Coordinate transformations.')
if name == 'utm':
return getUTMProjection(**kwargs)
elif name == 'RD83':
return pyproj.Proj(init="epsg:4745")
elif name == 'gk2':
return pyproj.Proj(init="epsg:31466")
elif name == 'gk3':
return pyproj.Proj(init="epsg:31467")
elif name == 'gk4':
return pyproj.Proj(init="epsg:31468")
elif name == 'gk5':
return pyproj.Proj(init="epsg:31469")
elif name == 'Soldner':
return pyproj.Proj(init="epsg:3068")
def _getXMLData(ele, name, default):
ret = default
if ele.getElementsByTagName(name):
ret = ele.getElementsByTagName(name)[0].childNodes[0].data
return ret
def _extractWPTS(wpts):
"""Handler for Waypoints in gpx xml-dom."""
w = dict()
for wpt in wpts:
if wpt.hasAttribute('lat'):
lat = float(wpt.getAttribute('lat'))
else:
continue
if wpt.hasAttribute('lon'):
lon = float(wpt.getAttribute('lon'))
else:
continue
name = _getXMLData(wpt, 'name', 'WP ' + str(len(w.keys())))
w[name] = {'lat': lat,
'lon': lon,
'time': _getXMLData(wpt, 'time', 0),
'desc': _getXMLData(wpt, 'desc', None),
}
return w
[docs]
def readGPX(fileName):
"""Extract GPS Waypoint from GPS Exchange Format (GPX).
Currently only simple waypoint extraction is supported.
<gpx version="1.0" creator="creator">
<metadata>
<name>Name</name>
</metadata>
<wpt lat="51." lon="11.">
<name>optional</name>
<time>optional</time>
<description>optional</description>
</wpt>
</gpx>
"""
from xml.dom.minidom import parse
dom = parse(fileName)
name = fileName
meta = dom.getElementsByTagName("metadata")
if len(meta) > 0:
name = _getXMLData(meta[0], 'name', fileName)
wpts = dom.getElementsByTagName("wpt")
return _extractWPTS(wpts), name
def readSimpleLatLon(filename, verbose=False):
"""Read Lat Lon coordinates from file.
Try converting automatic.
To be sure, provide format without "d" to ensure floating point format:
lon lat
or
marker lat lon
Returns
-------
list: []
lon lat name time
"""
def conv_(deg):
"""Convert degree into floating vpoint."""
ret = 0.0
if 'd' in deg:
# 10d???
d = deg.split('d')
ret = float(d[0])
if "'" in d[1]:
# 10d23'2323.44''
m = d[1].split("'")
if len(m) > 1:
ret += float(m[0]) / 60.
ret += float(m[1]) / 3600.
else:
# 10d23.2323
ret += float(d[1]) / 60.
else:
# 10.4432323
ret = float(deg)
return ret
# def conv_(...):
w = []
with open(filename, 'r') as fi:
content = fi.readlines()
for line in content:
if line[0] == '#':
continue
vals = line.split()
if len(vals) == 2: # lon lat
w.append((conv_(vals[1]), conv_(vals[0]), '', 'time'))
if len(vals) == 3:
# marker lat lon
w.append((conv_(vals[2]), conv_(vals[1]), vals[0], 'time'))
if verbose:
print(w[-1])
return w
def GK2toUTM(ea, no=None, zone=32):
"""Transform Gauss-Krueger zone 2 into UTM (for backward compatibility)."""
return GKtoUTM(ea, no, zone, gkzone=2)
def GK3toUTM(ea, no=None, zone=32):
"""Transform Gauss-Krueger zone 3 into UTM (for backward compatibility)."""
return GKtoUTM(ea, no, zone, gkzone=3)
def GK4toUTM(ea, no=None, zone=32):
"""Transform Gauss-Krueger zone 4 into UTM (for backward compatibility)."""
return GKtoUTM(ea, no, zone, gkzone=4)
def GK5toUTM(ea, no=None, zone=32):
"""Transform Gauss-Krueger zone 5 into UTM (for backward compatibility)."""
return GKtoUTM(ea, no, zone, gkzone=5)
[docs]
def GKtoUTM(ea, no=None, zone=32, gk=None, gkzone=None):
"""Transform any Gauss-Krueger to UTM autodetect GK zone from offset."""
if gk is None and gkzone is None:
if no is None:
rr = ea[0][0]
else:
if isinstance(ea, list) or isinstance(ea, tuple):
rr = ea[0]
else:
rr = ea
gkzone = int(floor(rr * 1e-6))
if gkzone <= 0 or gkzone >= 5:
print("cannot detect valid GK zone")
pyproj = pg.optImport('pyproj', 'coordinate transformations')
if pyproj is None:
return None
gk = pyproj.Proj(init="epsg:"+str(31464+gkzone))
wgs84 = pyproj.Proj(init="epsg:4326") # pure ellipsoid to double transform
utm = getUTMProjection(zone=zone, ellps='WGS84') # UTM
if no is None: # two-column matrix
lon, lat = pyproj.transform(gk, wgs84, ea[0], ea[1])
else:
lon, lat = pyproj.transform(gk, wgs84, ea, no)
return utm(lon, lat)
def convddmm(num):
"""Convert numeric position into degree and minute."""
dd = np.floor(num / 100.)
r1 = num - dd * 100.
return dd + r1 / 60.
def readGeoRefTIF(file_name):
"""Read geo-referenced TIFF file and return image and bbox.
plt.imshow(im, ext = bbox.ravel()), bbox might need transform.
"""
try:
import gdal
from osgeo.gdalconst import GA_ReadOnly
except ImportError:
sys.stderr.write("no module osgeo\n")
dataset = gdal.Open(file_name, GA_ReadOnly)
geotr = dataset.GetGeoTransform()
projection = dataset.GetProjection()
im = np.flipud(mpimg.imread(file_name))
tifx, tify, dx = geotr[0], geotr[3], geotr[1]
bbox = [[tifx, tifx + im.shape[1] * dx],
[tify - im.shape[0] * dx, tify]]
return im, bbox, projection
