UTM Conversion Utility

environment.yml

@@ -11,7 +11,7 @@ dependencies:
  - nose
  - numpy
  - numexpr
- - pyproj
+ - pyproj=2.2
  - psutil
  - obspy
  - vispy

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.017'
 
 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.run_geod(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

@@ -15,8 +15,10 @@
 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.017"
 
 
 class PH5Response(object):
@@ -68,6 +70,82 @@ def get_n_i(self, sensor_keys, datalogger_keys):
                 return ph5_resp.n_i
 
 
+class UTMConversions:
+
+    def __init__(self, lat, lon, side, zone):
+        self.epsg_wgs84 = "EPSG:4326"
+
+        if side is not None and lat is None:
+            if side == "S":
+                epsgroot = "327"
+            elif side == "N":
+                epsgroot = "326"
+        elif lat is not None:
+            if lat < 0.0:
+                epsgroot = "327"
+                side = 'S'
+            elif lat >= 0.0:
+                epsgroot = "326"
+                side = 'N'
+
+        if zone is not None and lon is None:
+            self.epsg_utm = "EPSG:" + epsgroot + str(zone)
+        elif lon is not None:
+            zone = self.lon2zone(lon)
+            self.epsg_utm = "EPSG:" + epsgroot + str(zone)
+
+        self.utm_zone = zone
+        self.hemisphere = side
+
+    def lon2zone(self, lon):
+        '''   Get UTM zone from longitude   '''
+        zone = int(math.floor((float(lon) + 180.)/6.0) % 60) + 1
+        return zone
+
+    def utm2latlong(self, easting, northing):
+        transformer = Transformer.from_crs(self.epsg_utm, self.epsg_wgs84,
+                                           always_xy=True)
+        lon, lat = transformer.transform(float(easting), float(northing))
+        return (lat, lon)
+
+    def geod2utm(self, lat, lon, elev):  # utility function
+        transformer = Transformer.from_crs(self.epsg_wgs84, self.epsg_utm,
+                                           always_xy=True)
+        easting, northing = transformer.transform(lon, lat)
+        return (northing, easting, elev)  # e.g. Y, X, Z
+
+    def latlong2utm(self, lat, lon):
+        transformer = Transformer.from_crs(self.epsg_wgs84, self.epsg_utm,
+                                           always_xy=True)
+        easting, northing = transformer.transform(lon, lat)
+        return (easting, northing, self.utm_zone, self.hemisphere)
+
+
+def tspc_lat_long(easting, northing):  # Texas State Plane Coords
+    epsg_wgs84 = "EPSG:4326"
+    epsg_sp4202 = "EPSG:2276"  # State Plane Coords in US FEET
+    transformer = Transformer.from_crs(epsg_sp4202, epsg_wgs84,
+                                       always_xy=True)
+    lon, lat = transformer.transform(easting, northing)
+    return (lon, lat)
+
+
+def run_geod(lat0, lon0, lat1, lon1, scalar=1.0):
+    ELLIPSOID = 'WGS84'
+
+    config = "+ellps={0}".format(ELLIPSOID)
+
+    g = Geod(config)
+
+    az, baz, dist = g.inv(lon0, lat0, lon1, lat1)
+
+    if dist:
+        dist /= scalar
+
+    # Return list containing azimuth, back azimuth, distance
+    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  # for geod2utm, run_geod
 from ph5.core import segy_h, ebcdic
 
-PROG_VERSION = '2018.268'
+PROG_VERSION = '2020.015'
 LOGGER = logging.getLogger(__name__)
 
 os.environ['TZ'] = 'UTC'
@@ -787,11 +786,11 @@ 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']
+                utmc = ph5utils.UTMConversions(lat, lon, None, None)
+                Y, X, Z = utmc.geod2utm(lat, lon, elev)
                 s, vx, vy = pick_values_32(X, Y)
 
                 tra['coordScale'] = s
@@ -844,11 +843,11 @@ 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']
+                    utmc = ph5utils.UTMConversions(lat, lon, None, None)
+                    Y, X, Z = utmc.geod2utm(lat, lon, elev)
 
                     s, vx, vy = pick_values_32(X, Y)
                     tra['sourceLongOrX'] = vx
@@ -873,10 +872,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.run_geod(lat1, lon1, lat2, lon2,
+                                                FACTS[UNITS])
                 tra['sourceToRecDist'] = az_baz_dist[2]
             except Exception as e:
                 # sys.stderr.write (e.message)
@@ -1058,23 +1061,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,12 @@
 from random import randint
 from inspect import stack
 from ast import literal_eval
-from pyproj import Geod
+from ph5.core import ph5utils  # for run_geod
 import simplekml as kml
 
 from ph5.core import timedoy
 
-PROG_VERSION = '2019.14'
+PROG_VERSION = '2020.015'
 LOGGER = logging.getLogger(__name__)
 
 FACTS = {'km': 1000., 'm': 1., 'dm': 1. / 10., 'cm': 1. / 100.,
@@ -216,15 +216,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.run_geod(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):
         '''