diff --git a/ph5/core/cs2cs.py b/ph5/core/cs2cs.py deleted file mode 100755 index d641792c..00000000 --- a/ph5/core/cs2cs.py +++ /dev/null @@ -1,92 +0,0 @@ -#!/usr/bin/env pnpython3 -# -# Convert geographic coordinates to UTM and back. -# Assumes reference ellipsoids are the same and are 3D. -# -# Steve Azevedo, January 2013 -# - -# https://code.google.com/p/pyproj/ -from pyproj import Proj - -# The proj4 program cs2cs must be in your path -# CS2CS = 'cs2cs' - - -def _sign(val, latlon): - ret = val - try: - nsew = str(val[0]) - nsew.upper() - if nsew == 'N' or nsew == 'S' or nsew == 'E' or nsew == 'W': - return ret - - if nsew == '+': - val = val[1:] - - if latlon == 'lat': - if nsew == '-': - ret = 'S' + val[1:] - else: - ret = 'N' + val - elif latlon == 'lon': - if nsew == '-': - ret = 'W' + val[1:] - else: - ret = 'E' + val - - except IndexError: - pass - - return ret - - -def lon2zone(lon): - ''' Get UTM zone from longitude, brute force method ''' - # zone 1 = -180 -> -174 - # zone 2 = -174 -> -168 - for zone in range(1, 60): - ebound = (zone * 6) - 180 - wbound = ebound - 6 - if lon >= wbound and lon <= ebound: - return zone - - return None - - -def utm2geod(zn, datum, X, Y, Z): - ''' Convert UTM coordinates to geodetic coordinates ''' - p = Proj(proj='utm', zone=zn, ellps=datum) - - lon, lan = p(X, Y, inverse=True) - - return lat, lon, Z - - -def geod2utm(zn, datum, lat, lon, elev): - ''' Convert geodetic coordinates to UTM ''' - if zn is None: - zn = lon2zone(lon) - - p = Proj(proj='utm', zone=zn, ellps=datum) - - X, Y = p(lon, lat) - - # Return Y, X, Z - return Y, X, elev - - -if __name__ == '__main__': - lat = 34.023786 - lon = -106.898492 - elev = 1456.0 - zone = 13 - - z = lon2zone(lon) - print z - - Y, X, Z = geod2utm(zone, 'WGS84', lat, lon, elev) - print Y, X, Z - - lat, lon, elev = utm2geod(zone, 'WGS84', X, Y, Z) - print lat, lon, elev diff --git a/ph5/core/ph5api.py b/ph5/core/ph5api.py index b71a310a..0f7d2c61 100755 --- a/ph5/core/ph5api.py +++ b/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.091' LOGGER = logging.getLogger(__name__) PH5VERSION = columns.PH5VERSION @@ -329,6 +328,7 @@ 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] @@ -336,11 +336,12 @@ def get_offset(self, sta_line, sta_id, evt_line, evt_id): 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.lat_lon_to_geod(lat0, lon0, + lat1, lon1, + FACTS_M[UNITS]) except Exception as e: LOGGER.warning("Couldn't get offset. {0}".format(repr(e))) return {} - return {'event_id_s': evt_id, 'receiver_id_s': sta_id, 'azimuth/value_f': az, 'azimuth/units_s': 'degrees', 'offset/value_d': dist, 'offset/units_s': 'm'} @@ -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 diff --git a/ph5/core/ph5utils.py b/ph5/core/ph5utils.py index c2c83f57..93e9247e 100644 --- a/ph5/core/ph5utils.py +++ b/ph5/core/ph5utils.py @@ -8,6 +8,7 @@ # 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 @@ -15,8 +16,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.091" class PH5Response(object): @@ -68,6 +71,105 @@ def get_n_i(self, sensor_keys, datalogger_keys): return ph5_resp.n_i +def lat_lon_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=False) + easting, northing = transformer.transform(lat, lon) + return (northing, easting, utm_zone, hemisphere) + + +def lat_lon_elev_to_utm(lat, lon, elev): + northing, easting, utm_zone, hemisphere = lat_lon_to_utm(lat, lon) + return (northing, easting, elev) + + +def utm_to_lat_lon(northing, easting, 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=False) + lat, lon = transformer.transform(float(easting), float(northing)) + return (lat, lon) + + +def tspc_lat_lon(northing, easting): + # 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=False) + lat, lon = transformer.transform(easting, northing) + return (lat, lon) + + +def lat_lon_to_geod(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 > 0.0: + dist /= scalar + return az, baz, dist + + +def lat_lon_to_ups_south(lat, lon): + # Universal Polar Stereographic, southern hemisphere + epsg_wgs84 = "EPSG:4326" + epsg_ups_south = "EPSG:32761" + transformer = Transformer.from_crs(epsg_wgs84, + epsg_ups_south, + always_xy=False) + northing, easting = transformer.transform(lat, lon) + return (northing, easting) + + +def lat_lon_to_ups_north(lat, lon): + # Universal Polar Stereographic, northern hemisphere + epsg_wgs84 = "EPSG:4326" + epsg_ups_north = "EPSG:32661" + transformer = Transformer.from_crs(epsg_wgs84, + epsg_ups_north, + always_xy=False) + northing, easting = transformer.transform(lat, lon) + return (northing, easting) + + +def ups_south_to_lat_lon(northing, easting): + # Universal Polar Stereogaphic, southern hemisphere + epsg_wgs84 = "EPSG:4326" + epsg_ups_south = "EPSG:32761" + transformer = Transformer.from_crs(epsg_ups_south, + epsg_wgs84, + always_xy=False) + lat, lon = transformer.transform(northing, easting) + return (lat, lon) + + +def ups_north_to_lat_lon(northing, easting): + # Universal Polar Stereogaphic, northern hemisphere + epsg_wgs84 = "EPSG:4326" + epsg_ups_north = "EPSG:32661" + transformer = Transformer.from_crs(epsg_ups_north, + epsg_wgs84, + always_xy=False) + lat, lon = transformer.transform(northing, easting) + return (lat, lon) + + """ =============== Utility methods diff --git a/ph5/core/segyfactory.py b/ph5/core/segyfactory.py index 066ad4fa..e51dfb52 100644 --- a/ph5/core/segyfactory.py +++ b/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.091' 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.lat_lon_elev_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.lat_lon_elev_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.lat_lon_to_geod(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 # diff --git a/ph5/entry_points.py b/ph5/entry_points.py index a059ebcf..80b539e6 100644 --- a/ph5/entry_points.py +++ b/ph5/entry_points.py @@ -5,7 +5,7 @@ # Dave Thomas, 2019-08-06 -PROG_VERSION = '2019.218' +PROG_VERSION = '2020.091' class EntryPointTypes(): @@ -277,5 +277,15 @@ 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), + EntryPoint('ups_latlong', + 'ph5.utilities.ups_latlong:main', + 'Converts UPS coordinates to/from' + ' Latitudes & Longitudes.', + type=EntryPointTypes.INGESTION), ] } diff --git a/ph5/utilities/novenqc.py b/ph5/utilities/novenqc.py index 64e6bcd4..658f71ed 100644 --- a/ph5/utilities/novenqc.py +++ b/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.091' + LOGGER = logging.getLogger(__name__) FACTS = {'km': 1000., 'm': 1., 'dm': 1. / 10., 'cm': 1. / 100., @@ -216,15 +217,14 @@ 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.lat_lon_to_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): ''' diff --git a/ph5/utilities/segd2ph5.py b/ph5/utilities/segd2ph5.py index f4856808..42df5278 100644 --- a/ph5/utilities/segd2ph5.py +++ b/ph5/utilities/segd2ph5.py @@ -15,11 +15,12 @@ import json import re from math import modf -from ph5.core import experiment, columns, segdreader +from ph5.core import experiment, columns, segdreader, ph5utils from ph5 import LOGGING_FORMAT -from pyproj import Proj, transform -PROG_VERSION = "2019.252" + +PROG_VERSION = "2020.085" + LOGGER = logging.getLogger(__name__) MAX_PH5_BYTES = 1073741824 * 100. # 100 GB (1024 X 1024 X 1024 X 2) @@ -922,50 +923,6 @@ def writeINDEX(): MAP_INFO = {} -def txncsptolatlon(northing, easting): - ''' - Sweetwater - Convert texas state plane coordinates in feet to - geographic coordinates, WGS84. - ''' - # Texas NC state plane feet Zone 4202 - sp = Proj(init='epsg:32038') - # WGS84, geographic - wgs = Proj(init='epsg:4326', proj='latlong') - # Texas SP coordinates: survey foot is 1200/3937 meters - lon, lat = transform(sp, wgs, easting * 0.30480060960121924, - northing * 0.30480060960121924) - - return lat, lon - - -def utmcsptolatlon(northing, easting): - ''' - Mount Saint Helens - Convert UTM to - geographic coordinates, WGS84. - ''' - # UTM - new_UTM = re.split(r'(\d+)', UTM) - utmzone = str(new_UTM[1]) - - if str(new_UTM[2]).upper() == 'N': - NS = 'north' - elif str(new_UTM[2]).upper() == 'S': - NS = 'south' - else: - NS = 'north' - - utmc = Proj("+proj=utm +zone="+utmzone+" +"+NS+" +ellps=WGS84") - print - # WGS84, geographic - wgs = Proj(init='epsg:4326', proj='latlong') - # - lon, lat = transform(utmc, wgs, easting, northing) - - return lat, lon - - def main(): import time then = time.time() @@ -1117,14 +1074,22 @@ def get_node(sd): try: if UTM: # UTM - LAT, LON = utmcsptolatlon( - SD.trace_headers.trace_header_N[ - 4].receiver_point_Y_final / 10., - SD.trace_headers.trace_header_N[ - 4].receiver_point_X_final / 10.) + new_UTM = re.split(r'(\d+)', UTM) + utmzone = str(new_UTM[1]) + + NS = str(new_UTM[2]).upper() + if NS != 'N' and NS != 'S': + NS = 'N' + easting = SD.trace_headers.trace_header_N[ + 4].receiver_point_X_final / 10. + northing = SD.trace_headers.trace_header_N[ + 4].receiver_point_Y_final / 10. + LAT, LON = ph5utils.utm_to_lat_lon( + northing, easting, + NS, utmzone) elif TSPF: # Texas State Plane coordinates - LAT, LON = txncsptolatlon( + LAT, LON = ph5utils.tspc_lat_lon( SD.trace_headers.trace_header_N[ 4].receiver_point_Y_final / 10., SD.trace_headers.trace_header_N[ diff --git a/ph5/utilities/tests/test_utmtolatlong.py b/ph5/utilities/tests/test_utmtolatlong.py new file mode 100644 index 00000000..6471f0a2 --- /dev/null +++ b/ph5/utilities/tests/test_utmtolatlong.py @@ -0,0 +1,104 @@ +import unittest +import pyproj +from ph5.core.ph5utils import lat_lon_to_utm +from ph5.core.ph5utils import utm_to_lat_lon +from ph5.core.ph5utils import lat_lon_elev_to_utm +from ph5.core.ph5utils import tspc_lat_lon +from ph5.core.ph5utils import lat_lon_to_geod +from ph5.core.ph5utils import lat_lon_to_ups_north +from ph5.core.ph5utils import lat_lon_to_ups_south +from ph5.core.ph5utils import ups_north_to_lat_lon +from ph5.core.ph5utils import ups_south_to_lat_lon + + +PROG_VERSION = '2020.091' + + +class TestUTMconversion(unittest.TestCase): + + def test_is_valid_utm_conversion_north(self): + # location: PIC + lat, lon = utm_to_lat_lon(3771967, 322916, 'N', 13) + self.assertAlmostEqual(lat, 34.07349577107704) + self.assertAlmostEqual(lon, -106.91909595147378) + + def test_is_valid_utm_conversion_south(self): + # location: Castle Rock, Antarctica + lat, lon = utm_to_lat_lon(1361594, 540947, 'S', 58) + self.assertAlmostEqual(lat, -77.81567398301094) + self.assertAlmostEqual(lon, 166.73816638798527) + + def test_for_correct_type(self): + # heck, no are invalid eastings/northings + with self.assertRaises(ValueError): + lat, lon = utm_to_lat_lon('heck', 'no', 'S', 58) + + def test_for_correct_value(self): + # 666 is invalid UTM zone + with self.assertRaises(pyproj.exceptions.CRSError): + lat, lon = utm_to_lat_lon(1361594, 540947, 'S', 666) + + def test_is_valid_geod_conversion(self): + # location: PIC + northing, easting, elev =\ + lat_lon_elev_to_utm(34.0734958, -106.9190959, 1456.0) + self.assertAlmostEqual(northing, 3771967.003118457) + self.assertAlmostEqual(easting, 322916.0048106084) + self.assertAlmostEqual(elev, 1456.0) + + def test_is_valid_utm_latlong_conversion_north(self): + # location: PIC + northing, easting, zone, hemisphere =\ + lat_lon_to_utm(34.0734958, -106.9190959) + self.assertAlmostEqual(northing, 3771967.003118457) + self.assertAlmostEqual(easting, 322916.0048106084) + self.assertEqual(zone, 13) + self.assertAlmostEqual(hemisphere, 'N') + + def test_is_valid_utm_latlong_conversion_south(self): + # location: Castle Rock Antarctica + northing, easting, zone, hemisphere =\ + lat_lon_to_utm(-77.81567398301094, 166.73816638798527) + self.assertAlmostEqual(northing, 1361594.0) + self.assertAlmostEqual(easting, 540947.0) + self.assertEqual(zone, 58) + self.assertAlmostEqual(hemisphere, 'S') + + def test_is_valid_tsp_conversion(self): + # Sweetwater, Texas, units US FEET + lat, lon = tspc_lat_lon(6858888, 1380811) + self.assertAlmostEqual(lat, 32.468972700219375) + self.assertAlmostEqual(lon, -100.40568335900281) + + def test_is_valid_lat_lon_to_geod_conversion(self): + # socorro to albuquerque in miles + az, baz, dist = lat_lon_to_geod(34.0543, -106.907, + 35.1053, -106.646, + 1609.344) + self.assertAlmostEqual(az, 11.533074930503515) + self.assertAlmostEqual(baz, -168.3187871800244) + self.assertAlmostEqual(dist, 73.95826059337135) + + def test_is_valid_ups_conversion_north(self): + lat, lon = ups_north_to_lat_lon(1500000.0, 2500000.0) + self.assertAlmostEqual(lat, 83.63731756105707) + self.assertAlmostEqual(lon, 45.0) + + def test_is_valid_ups_conversion_south(self): + lat, lon = ups_south_to_lat_lon(1500000.0, 2500000.0) + self.assertAlmostEqual(lat, -83.63731756105707) + self.assertAlmostEqual(lon, 135.0) + + def test_is_valid_ups_latlong_conversion_north(self): + northing, easting = lat_lon_to_ups_north(83.63731756105707, 135.) + self.assertAlmostEqual(northing, 2500000.0) + self.assertAlmostEqual(easting, 2500000.0) + + def test_is_valid_ups_latlong_conversion_south(self): + northing, easting = lat_lon_to_ups_south(-83.63731756105707, 45.) + self.assertAlmostEqual(northing, 2500000.0) + self.assertAlmostEqual(easting, 2500000.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/ph5/utilities/ups_latlong.py b/ph5/utilities/ups_latlong.py new file mode 100644 index 00000000..e2adbaf4 --- /dev/null +++ b/ph5/utilities/ups_latlong.py @@ -0,0 +1,184 @@ +# This script is used to convert UPS coordinates to/from Latitude/Longitude +# Dave Thomas, 2020-03-27 + +from __future__ import (print_function) +from ph5.core import ph5utils +import argparse +import logging + +PROG_VERSION = '2020.091' +LOGGER = logging.getLogger(__name__) + + +def convert_file_from_ups(infile, outfile): + # batch file method + counter = 0 + with open(outfile, "w") as theoutfile, open(infile) as theinfile: + ws = "Hemisphere, Northing, Easting, => Latitude, Longitude\n" + theoutfile.write(ws) + for line in theinfile: + counter = counter + 1 + if counter == 1: + # skip header line + continue + fieldset = line.strip().split(",") + if len(fieldset) <= 1: + continue + hemisphere = str(fieldset[0]).upper() + northing = float(fieldset[1]) + easting = float(fieldset[2]) + try: + if hemisphere == 'N': + lat, lon = ph5utils.ups_north_to_lat_lon(northing, easting) + else: + # default to southern hemisphere + lat, lon = ph5utils.ups_south_to_lat_lon(northing, easting) + ps = "hemisphere=%s, northing=%.1f, easting=%.1f" + printstr1 = (ps % (hemisphere, northing, easting)) + printstr2 = "Lat = %.7f, Lon = %.7f" % (lat, lon) + LOGGER.info(printstr1 + " => " + printstr2) + outstr = "%s, %.1f, %.1f, %.7f, %.7f" \ + % (hemisphere, northing, easting, lat, lon) + theoutfile.write(outstr + "\n") + except ValueError: + LOGGER.error("There was an error in ups data conversion.") + return + + +def convert_file_from_latlong(infile, outfile): + # batch file method + counter = 0 + with open(outfile, "w") as theoutfile, open(infile) as theinfile: + ws = "Latitude, Longitude, => Hemisphere, Northing, Easting\n" + theoutfile.write(ws) + for line in theinfile: + counter = counter + 1 + if counter == 1: + # skip header line + continue + fieldset = line.strip().split(",") + if len(fieldset) <= 1: + continue + lat = float(fieldset[0]) + lon = float(fieldset[1]) + try: + if lat > 0.0: + hemisphere = 'N' + northing, easting =\ + ph5utils.lat_lon_to_ups_north(lat, lon) + else: + hemisphere = 'S' + northing, easting =\ + ph5utils.lat_lon_to_ups_south(lat, lon) + ps = "hemisphere=%s, northing=%.1f, easting=%.1f" + printstr1 = "Lat = %.7f, Lon = %.7f" % (lat, lon) + printstr2 = (ps % (hemisphere, northing, easting)) + LOGGER.info(printstr1 + " => " + printstr2) + outstr = "%.7f, %.7f, %s, %.1f, %.1f" \ + % (lat, lon, hemisphere, northing, easting) + theoutfile.write(outstr + "\n") + except ValueError: + LOGGER.error("There was an error in ups data conversion.") + return + + +def doinline_from_ups(easting, northing, hemisphere='S'): + # interactive method + try: + indat = ("Input: easting=%.1f, northing=%.1f, hemisphere=%s" + % (float(easting), float(northing), hemisphere)) + LOGGER.info(indat) + + if hemisphere == 'N': + lat, lon = ph5utils.ups_north_to_lat_lon(northing, easting) + else: + lat, lon = ph5utils.ups_south_to_lat_lon(northing, easting) + outstr = "Lat = %.7f, Lon = %.7f" % (lat, lon) + LOGGER.info(outstr) + LOGGER.info(''.join('*'*50)) + except ValueError: + LOGGER.error("There was an error in ups data conversion.") + return + + +def doinline_from_latlong(lat, lon): + # interactive method + try: + if abs(float(lat)) > 90.0: + LOGGER.error( + "Your latitude, {0}, is out of limit." + .format(lat)) + return + if float(lat) > 0.0: + hemisphere = 'Northern' + northing, easting = \ + ph5utils.lat_lon_to_ups_north(float(lat), float(lon)) + else: + hemisphere = 'Southern' + northing, easting = \ + ph5utils.lat_lon_to_ups_south(float(lat), float(lon)) + indat = "Input: lat=%.7f, lon=%.7f" % (float(lat), float(lon)) + LOGGER.info(indat) + outstr = "%s hemisphere, northing =%.1f, easting =%.1f" \ + % (hemisphere, float(northing), float(easting)) + LOGGER.info(outstr) + LOGGER.info(''.join('*'*50)) + except ValueError: + LOGGER.error("There was an error in ups data conversion.") + return + + +def main(): + desc = 'Convert UPS locations to/from longitudes, latitudes, \ +in batches (-i, -o) or inline (-e,-n,-s), (-x,-y)' + parser = argparse.ArgumentParser(description=desc) + group1 = parser.add_mutually_exclusive_group() + group1.add_argument("-u", "--ups", help="Convert UPS to Lat/Long", + action="store_true") + group1.add_argument("-l", "--latlong", help="Convert Lat/Long to UPS", + action="store_true") + parser.add_argument("-i", "--infile", help="Input file with UPS \ +Easting, Northing, Hemisphere, or Lats and Longs") + parser.add_argument("-o", "--outfile", help="Output file with \ +converted longitudes, latitudes or ups Coordinates") + parser.add_argument("-e", "--easting", help="Enter UPS Easting") + parser.add_argument("-n", "--northing", help="Enter UPS Northing") + parser.add_argument("-s", "--side", help="Side=Hemisphere, N or S") + parser.add_argument("-x", "--longitude", help="Enter longitude") + parser.add_argument("-y", "--latitude", help="Enter latitude") + + args = parser.parse_args() + + # for batch file method + if args.infile is not None and args.outfile is not None: + if args.ups == 1: + convert_file_from_ups(args.infile, args.outfile) + elif args.latlong == 1: + convert_file_from_latlong(args.infile, args.outfile) + else: + LOGGER.error("You must specify either -u/--ups, OR -l/--latlong") + + # for direct method + else: + if args.ups == 1: + if args.easting is not None and args.northing is not None: + if args.side is not None: + doinline_from_ups(args.easting, args.northing, + args.side.upper()) + else: + # default = Southern Hemisphere + doinline_from_ups(args.easting, args.northing, + 'S') + else: + LOGGER.error("You must specify northing and easting") + elif args.latlong == 1: + if args.longitude is not None and args.latitude is not None: + doinline_from_latlong(args.latitude, args.longitude) + else: + LOGGER.error("You must specify longitude and latitude.") + else: + LOGGER.error("You must choose -u/--ups, OR -l/--latlong") + + +if __name__ == '__main__': + main() diff --git a/ph5/utilities/utm_latlong.py b/ph5/utilities/utm_latlong.py new file mode 100644 index 00000000..24be35dd --- /dev/null +++ b/ph5/utilities/utm_latlong.py @@ -0,0 +1,165 @@ +# This script is used to convert UTM coordinates to/from Latitude/Longitude +# Dave Thomas, 2020-02-24 + +from __future__ import (print_function) +from ph5.core import ph5utils +import argparse +import logging + +PROG_VERSION = '2020.091' +LOGGER = logging.getLogger(__name__) + + +def convert_file_from_utm(infile, outfile): + # batch file method + counter = 0 + with open(outfile, "w") as theoutfile, open(infile) as theinfile: + ws = "Northing, Easting, Zone, Hemisphere, => Latitude, Longitude\n" + theoutfile.write(ws) + for line in theinfile: + counter = counter + 1 + if counter == 1: + # skip header line + continue + fieldset = line.strip().split(",") + if len(fieldset) <= 1: + continue + easting = float(fieldset[0]) + northing = float(fieldset[1]) + zone = int(fieldset[2]) + hemisphere = str(fieldset[3]).upper() + try: + lat, lon = ph5utils.utm_to_lat_lon(northing, easting, + hemisphere, zone) + ps = "northing=%.1f, easting=%.1f, zone=%d, hemisphere=%s" + printstr1 = (ps % (northing, easting, zone, hemisphere)) + printstr2 = "Lat = %.7f, Lon = %.7f" % (lat, lon) + LOGGER.info(printstr1 + " => " + printstr2) + outstr = "%.1f, %.1f, %d, %s, %.7f, %.7f" \ + % (northing, easting, zone, hemisphere, lat, lon) + theoutfile.write(outstr + "\n") + except ValueError: + LOGGER.error("There was an error in UTM data conversion.") + return + + +def convert_file_from_latlong(infile, outfile): + # batch file method + counter = 0 + with open(outfile, "w") as theoutfile, open(infile) as theinfile: + ws = "Latitude, Longitude, =>Northing, Easting, Zone, Hemisphere\n" + theoutfile.write(ws) + for line in theinfile: + counter = counter + 1 + if counter == 1: + # skip header line + continue + fieldset = line.strip().split(",") + if len(fieldset) <= 1: + continue + lat = float(fieldset[0]) + lon = float(fieldset[1]) + try: + northing, easting, zone, hemisphere =\ + ph5utils.lat_lon_to_utm(lat, lon) + ps = "northing=%.1f, easting=%.1f, zone=%d, hemisphere=%s" + printstr1 = "Lat = %.7f, Lon = %.7f" % (lat, lon) + printstr2 = (ps % (northing, easting, zone, hemisphere)) + LOGGER.info(printstr1 + " => " + printstr2) + outstr = "%.7f, %.7f, %.1f, %.1f, %d, %s" \ + % (lat, lon, northing, easting, zone, hemisphere) + theoutfile.write(outstr + "\n") + except ValueError: + LOGGER.error("There was an error in UTM data conversion.") + return + + +def doinline_from_utm(easting, northing, zone, hemisphere): + # interactive method + try: + indat = ("Input: easting=%.1f, northing=%.1f, zone=%d, hemisphere=%s" + % (float(easting), float(northing), int(zone), hemisphere)) + LOGGER.info(indat) + + lat, lon = ph5utils.utm_to_lat_lon(northing, easting, + hemisphere, zone) + outstr = "Lat = %.7f, Lon = %.7f" % (lat, lon) + LOGGER.info(outstr) + LOGGER.info(''.join('*'*50)) + except ValueError: + LOGGER.error("There was an error in UTM data conversion.") + return + + +def doinline_from_latlong(lat, lon): + # interactive method + try: + if abs(float(lat)) > 90.0: + LOGGER.error( + "Your latitude, {0}, is out of limit." + .format(lat)) + return + northing, easting, zone, hemisphere = \ + ph5utils.lat_lon_to_utm(float(lat), float(lon)) + indat = "Input: lat=%.7f, lon=%.7f" % (float(lat), float(lon)) + LOGGER.info(indat) + outstr = "northing =%.1f, easting =%.1f, zone =%s, hemisphere =%s" \ + % (float(northing), float(easting), zone, hemisphere) + LOGGER.info(outstr) + LOGGER.info(''.join('*'*50)) + except ValueError: + LOGGER.error("There was an error in UTM data conversion.") + return + + +def main(): + desc = 'Convert UTM locations to/from longitudes, latitudes, \ +in batches (-i, -o) or inline (-e,-n,-z,-s), (-x,-y)' + parser = argparse.ArgumentParser(description=desc) + group1 = parser.add_mutually_exclusive_group() + group1.add_argument("-u", "--utm", help="Convert UTM to Lat/Long", + action="store_true") + group1.add_argument("-l", "--latlong", help="Convert Lat/Long to UTM", + action="store_true") + parser.add_argument("-i", "--infile", help="Input file with UTM \ +Easting, Northing, Zone, Hemisphere, or Lats and Longs") + parser.add_argument("-o", "--outfile", help="Output file with \ +converted longitudes, latitudes or UTM Coordinates") + parser.add_argument("-e", "--easting", help="Enter UTM Easting") + parser.add_argument("-n", "--northing", help="Enter UTM Northing") + parser.add_argument("-z", "--zone", help="Enter UTM Zone") + parser.add_argument("-s", "--side", help="Side=Hemisphere, N or S") + parser.add_argument("-x", "--longitude", help="Enter longitude") + parser.add_argument("-y", "--latitude", help="Enter latitude") + + args = parser.parse_args() + + # for batch file method + if args.infile is not None and args.outfile is not None: + if args.utm == 1: + convert_file_from_utm(args.infile, args.outfile) + elif args.latlong == 1: + convert_file_from_latlong(args.infile, args.outfile) + else: + LOGGER.error("You must specify either -u/--utm, OR -l/--latlong") + + # for direct method + else: + if args.utm == 1: + if args.easting is not None and args.northing is not None \ + and args.zone is not None and args.side is not None: + doinline_from_utm(args.easting, args.northing, + args.zone, args.side.upper()) + else: + LOGGER.error("You must specify northing,easting,zone, & side") + elif args.latlong == 1: + if args.longitude is not None and args.latitude is not None: + doinline_from_latlong(args.latitude, args.longitude) + else: + LOGGER.error("You must specify longitude and latitude.") + else: + LOGGER.error("You must choose -u/--utm, OR -l/--latlong") + + +if __name__ == '__main__': + main() diff --git a/setup.py b/setup.py index 393efe2e..19404a75 100644 --- a/setup.py +++ b/setup.py @@ -54,7 +54,7 @@ 'nose', 'numpy', 'numexpr', - 'pyproj', + 'pyproj>=2', 'psutil', 'obspy', 'vispy',