UTM Conversion Utility

ph5/core/ph5api.py

@@ -11,11 +11,10 @@
 import re
 import numpy as np
 import math
-from pyproj import Geod
-from ph5.core import columns, experiment, timedoy
+from ph5.core import columns, experiment, timedoy, ph5utils
 from tables.exceptions import NoSuchNodeError
 
-PROG_VERSION = '2019.93'
+PROG_VERSION = '2020.066'
 
 LOGGER = logging.getLogger(__name__)
 PH5VERSION = columns.PH5VERSION
@@ -329,14 +328,16 @@ def get_offset(self, sta_line, sta_id, evt_line, evt_id):
             LOGGER.warning("Couldn't get offset.")
             return {}
         try:
+            UNITS = 'm'
             if sta_line in self.Array_t and evt_line in self.Event_t:
                 array_t = self.Array_t[sta_line]['byid'][sta_id][c]
                 event_t = self.Event_t[evt_line]['byid'][evt_id]
                 lon0 = array_t[0]['location/X/value_d']
                 lat0 = array_t[0]['location/Y/value_d']
                 lon1 = event_t['location/X/value_d']
                 lat1 = event_t['location/Y/value_d']
-                az, baz, dist = run_geod(lat0, lon0, lat1, lon1)
+                az, baz, dist = ph5utils.latlon2geod(lat0, lon0, lat1, lon1,
+                                                     FACTS_M[UNITS])
         except Exception as e:
             LOGGER.warning("Couldn't get offset. {0}".format(repr(e)))
             return {}
@@ -1664,23 +1665,6 @@ def by_id(rows, key='id_s', secondary_key=None, unique_key=True):
     return byid, order
 
 
-def run_geod(lat0, lon0, lat1, lon1):
-    UNITS = 'm'
-    ELLIPSOID = 'WGS84'
-
-    config = "+ellps={0}".format(ELLIPSOID)
-
-    g = Geod(config)
-
-    az, baz, dist = g.inv(lon0, lat0, lon1, lat1)
-
-    if dist:
-        dist /= FACTS_M[UNITS]
-
-    # Return list containing azimuth, back azimuth, distance
-    return az, baz, dist
-
-
 def rect(r, w, deg=0):
     # Convert from polar to rectangular coordinates
     # radian if deg=0; degree if deg=1

ph5/core/ph5utils.py

@@ -8,15 +8,18 @@
 # version: 2019.032
 # author: Lan Dam
 # add inday_breakup() function
+# Dave Thomas added UTM / Lat Long routines,consolidated all pyproj here
 
 import fnmatch
 from datetime import datetime, timedelta
 from obspy.geodetics import locations2degrees
 from ph5.core.timedoy import epoch2passcal, passcal2epoch, TimeDOY, TimeError
 import time
 import re
+from pyproj import Transformer, Geod
+import math
 
-PROG_VERSION = "2019.81"
+PROG_VERSION = "2020.066"
 
 
 class PH5Response(object):
@@ -68,6 +71,76 @@ def get_n_i(self, sensor_keys, datalogger_keys):
                 return ph5_resp.n_i
 
 
+def lat_long_to_utm(lat, lon):
+
+    epsg_wgs84 = "EPSG:4326"
+
+    if lat < 0.0:
+        epsgroot = "327"
+        hemisphere = 'S'
+    elif lat >= 0.0:
+        epsgroot = "326"
+        hemisphere = 'N'
+
+    utm_zone = int(math.floor((float(lon) + 180.)/6.0) % 60) + 1
+    epsg_utm = "EPSG:" + epsgroot + str(utm_zone)
+    transformer = Transformer.from_crs(epsg_wgs84,
+                                       epsg_utm,
+                                       always_xy=True)
+    easting, northing = transformer.transform(lon, lat)
+
+    return (easting, northing, utm_zone, hemisphere)
+
+
+def geod_to_utm(lat, lon, elev):
+
+    easting, northing, utm_zone, hemisphere = lat_long_to_utm(lat, lon)
+
+    return (northing, easting, elev)
+
+
+def utm_to_lat_long(easting, northing, hemisphere, zone):
+
+    epsg_wgs84 = "EPSG:4326"
+
+    if hemisphere == "S":
+        epsgroot = "327"
+    elif hemisphere == "N":
+        epsgroot = "326"
+
+    epsg_utm = "EPSG:" + epsgroot + str(zone)
+
+    transformer = Transformer.from_crs(epsg_utm, epsg_wgs84, always_xy=True)
+    lon, lat = transformer.transform(float(easting), float(northing))
+    return (lat, lon)
+
+
+def tspc_lat_long(easting, northing):
+    # Texas State Plane Coords, in US FEET
+    epsg_wgs84 = "EPSG:4326"
+    epsg_sp4202 = "EPSG:2276"
+    transformer = Transformer.from_crs(epsg_sp4202, epsg_wgs84,
+                                       always_xy=True)
+    lon, lat = transformer.transform(easting, northing)
+    return (lon, lat)
+
+
+def latlon2geod(lat0, lon0, lat1, lon1, scalar=1.0):
+    # Return azimuth, back azimuth, geodetic distance
+    ELLIPSOID = 'WGS84'
+
+    config = "+ellps={0}".format(ELLIPSOID)
+
+    g = Geod(config)
+
+    az, baz, dist = g.inv(lon0, lat0, lon1, lat1)
+
+    if dist:
+        dist /= scalar
+
+    return az, baz, dist
+
+
 """
 ===============
 Utility methods

ph5/core/segyfactory.py

@@ -16,11 +16,10 @@
 import string
 import sys
 import logging
-from pyproj import Geod
-from ph5.core.cs2cs import geod2utm
+from ph5.core import ph5utils
 from ph5.core import segy_h, ebcdic
 
-PROG_VERSION = '2018.268'
+PROG_VERSION = '2020.066'
 LOGGER = logging.getLogger(__name__)
 
 os.environ['TZ'] = 'UTC'
@@ -787,11 +786,10 @@ def set_trace_header(self):
 
         if self.utm is True:
             try:
-                Y, X, Z = geod2utm(None,  # Zone goes here
-                                   "WGS84",
-                                   self.array_t['location/Y/value_d'],
-                                   self.array_t['location/X/value_d'],
-                                   self.array_t['location/Z/value_d'])
+                lat = self.array_t['location/Y/value_d']
+                lon = self.array_t['location/X/value_d']
+                elev = self.array_t['location/Z/value_d']
+                Y, X, Z = ph5utils.geod_to_utm(lat, lon, elev)
                 s, vx, vy = pick_values_32(X, Y)
 
                 tra['coordScale'] = s
@@ -844,11 +842,10 @@ def set_trace_header(self):
 
             if self.utm:
                 try:
-                    Y, X, Z = geod2utm(None,  # Zone goes here
-                                       "WGS84",
-                                       self.event_t['location/Y/value_d'],
-                                       self.event_t['location/X/value_d'],
-                                       self.event_t['location/Z/value_d'])
+                    lat = self.event_t['location/Y/value_d']
+                    lon = self.event_t['location/X/value_d']
+                    elev = self.event_t['location/Z/value_d']
+                    Y, X, Z = ph5utils.geod_to_utm(lat, lon, elev)
 
                     s, vx, vy = pick_values_32(X, Y)
                     tra['sourceLongOrX'] = vx
@@ -873,10 +870,14 @@ def set_trace_header(self):
             tra['sourceToRecDist'] = self.offset_t['offset/value_d']
         else:
             try:
-                az_baz_dist = run_geod(self.event_t['location/Y/value_d'],
-                                       self.event_t['location/X/value_d'],
-                                       self.array_t['location/Y/value_d'],
-                                       self.array_t['location/X/value_d'])
+                lat1 = self.event_t['location/Y/value_d']
+                lon1 = self.event_t['location/X/value_d']
+                lat2 = self.array_t['location/Y/value_d']
+                lon2 = self.array_t['location/X/value_d']
+
+                UNITS = 'm'
+                az_baz_dist = ph5utils.latlon2geod(lat1, lon1, lat2, lon2,
+                                                   FACTS[UNITS])
                 tra['sourceToRecDist'] = az_baz_dist[2]
             except Exception as e:
                 # sys.stderr.write (e.message)
@@ -1058,23 +1059,6 @@ def pick_values_32(X, Y):
     return sx, vx, vy
 
 
-def run_geod(lat0, lon0, lat1, lon1):
-    ELLIPSOID = 'WGS84'
-    UNITS = 'm'
-
-    config = "+ellps={0}".format(ELLIPSOID)
-
-    g = Geod(config)
-
-    az, baz, dist = g.inv(lon0, lat0, lon1, lat1)
-
-    if dist:
-        dist /= FACTS[UNITS]
-
-    # Return list containing azimuth, back azimuth, distance
-    return az, baz, dist
-
-
 #
 # Write standard SEG-Y reel header
 #

ph5/entry_points.py

@@ -5,7 +5,7 @@
 
 # Dave Thomas, 2019-08-06
 
-PROG_VERSION = '2019.218'
+PROG_VERSION = '2019.330'
 
 
 class EntryPointTypes():
@@ -277,5 +277,10 @@ def __init__(self):
                            'Takes PH5 files and returns time series '
                            'availability info.',
                            type=EntryPointTypes.CLIENT),
+                EntryPoint('utm_latlong',
+                           'ph5.utilities.utm_latlong:main',
+                           'Converts UTM coordinates to/from'
+                           ' Latitudes & Longitudes.',
+                           type=EntryPointTypes.INGESTION),
                 ]
             }

ph5/utilities/novenqc.py

@@ -12,12 +12,13 @@
 from random import randint
 from inspect import stack
 from ast import literal_eval
-from pyproj import Geod
+from ph5.core import ph5utils
 import simplekml as kml
 
 from ph5.core import timedoy
 
-PROG_VERSION = '2019.14'
+PROG_VERSION = '2020.066'
+
 LOGGER = logging.getLogger(__name__)
 
 FACTS = {'km': 1000., 'm': 1., 'dm': 1. / 10., 'cm': 1. / 100.,
@@ -216,15 +217,13 @@ def qc_dist(ys, xs, zs):
         # Need to check for UTM and convert to lat/lon
         #
         units = 'm'
-        ellipsoid = 'WGS84'
-        config = "+ellps={0}".format(ellipsoid)
-        g = Geod(config)
-        az, baz, dist = g.inv(xs[0], ys[0], xs[1], ys[1])
+        az, baz, dist = ph5utils.latlon2geod(ys[0], xs[0], ys[1], xs[1],
+                                             FACTS[units])
         if len(zs) > 1:
             zdelta = float(zs[1]) - float(zs[0])
         else:
             zdelta = 0.0
-        return az, baz, dist / FACTS[units], zdelta
+        return az, baz, dist, zdelta
 
     def qc_time(t1, t2):
         '''