From c3bb1a3a9d51015e08c3efa2b8104595b1fdafc7 Mon Sep 17 00:00:00 2001 From: shaunwbell Date: Tue, 6 Feb 2024 15:55:23 -0800 Subject: [PATCH] Create find_closest_ctd.py (#297) * Create find_closest_ctd.py * add haversine * Update find_closest_ctd.py * Update find_closest_ctd.py * more cleanup * Update find_closest_ctd.py --- src/ecofocipy/math/haversine.py | 45 ++++++ tools/find_closest_ctd.py | 252 ++++++++++++++++++++++++++++++++ 2 files changed, 297 insertions(+) create mode 100644 src/ecofocipy/math/haversine.py create mode 100644 tools/find_closest_ctd.py diff --git a/src/ecofocipy/math/haversine.py b/src/ecofocipy/math/haversine.py new file mode 100644 index 0000000..794711c --- /dev/null +++ b/src/ecofocipy/math/haversine.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python +# Haversine formula example in Python +# Author: Wayne Dyck + +import math + +import numpy as np + + +def distance(origin, destination): + lat1, lon1 = origin + lat2, lon2 = destination + radius = 6371 # km + + dlat = math.radians(lat2-lat1) + dlon = math.radians(lon2-lon1) + a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \ + * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2) + c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)) + d = radius * c + + return d + +def nearest_point(origin, latpoints, lonpoints, grid='1d'): + + if grid == '1d': + dist = np.zeros((np.shape(latpoints)[0], np.shape(lonpoints)[0])) + + for i, lat in enumerate(latpoints): + for j, lon in enumerate(lonpoints): + dist[i,j] = distance(origin,[lat,lon]) + + lati, loni = np.where(dist == dist.min()) + return (dist.min(), latpoints[lati[0]], lonpoints[loni[0]], lati[0], loni[0] ) + + elif grid == '2d': + dist = np.zeros_like(latpoints) + + for i, latrow in enumerate(latpoints): + for j, lonrow in enumerate(latrow): + dist[i,j] = distance(origin,[latpoints[i,j],lonpoints[i,j]]) + + lati, loni = np.where(dist == dist.min()) + + return (dist.min(), latpoints[lati[0]][loni[0]], lonpoints[lati[0]][loni[0]], lati[0], loni[0] ) \ No newline at end of file diff --git a/tools/find_closest_ctd.py b/tools/find_closest_ctd.py new file mode 100644 index 0000000..ab88b11 --- /dev/null +++ b/tools/find_closest_ctd.py @@ -0,0 +1,252 @@ +#!/usr/bin/env python + +""" + FindClosestCTD.py + + Using the EcoFOCI ruise log database, search for CTD's within a specified range of a location + + + History + ======= + 2024-02-06: Migrate to EcoFOCIpy Tools + 2019-07-15: Make python3 compliant: WIP + + Future + ====== + + Compatibility: + ============== + python >=3.8 **Tested + + """ + +import argparse +import datetime +import sys + +import mysql.connector +import numpy as np +from ecofocipy.math.haversine import sphered + +__author__ = "Shaun Bell" +__email__ = "shaun.bell@noaa.gov" +__created__ = datetime.datetime(2016, 9, 28) +__modified__ = datetime.datetime(2024, 2, 6) +__version__ = "0.1.1" +__status__ = "Development" + + +"""--------------------------------SQL Init----------------------------------------""" + + +class NumpyMySQLConverter(mysql.connector.conversion.MySQLConverter): + """ A mysql.connector Converter that handles Numpy types """ + + def _float32_to_mysql(self, value): + if np.isnan(value): + return None + return float(value) + + def _float64_to_mysql(self, value): + if np.isnan(value): + return None + return float(value) + + def _int32_to_mysql(self, value): + if np.isnan(value): + return None + return int(value) + + def _int64_to_mysql(self, value): + if np.isnan(value): + return None + return int(value) + + +def connect_to_DB(**kwargs): + # Open database connection + try: + db = mysql.connector.connect(use_pure=True, **kwargs) + except mysql.connector.Error as err: + """ + if err.errno == errorcode.ER_ACCESS_DENIED_ERROR: + print("Something is wrong with your user name or password") + elif err.errno == errorcode.ER_BAD_DB_ERROR: + print("Database does not exist") + else: + print(err) + """ + print("error - will robinson") + + db.set_converter_class(NumpyMySQLConverter) + + # prepare a cursor object using cursor() method + cursor = db.cursor(dictionary=True) + prepcursor = db.cursor(prepared=True) + return (db, cursor) + + +def close_DB(db): + # disconnect from server + db.close() + + +def read_data(db, cursor, table, yearrange): + + """ + + """ + sql = ( + "SELECT `id`,`LatitudeDeg`,`LatitudeMin`,`LongitudeDeg`," + "`LongitudeMin`,`ConsecutiveCastNo`,`UniqueCruiseID`,`GMTDay`," + "`GMTMonth`,`GMTYear`,`MaxDepth` from `{table}` WHERE `GMTYear` BETWEEN '{startyear}' AND '{endyear}'" + ).format(table=table, startyear=yearrange[0], endyear=yearrange[1]) + print(sql) + + result_dic = {} + try: + # Execute the SQL command + cursor.execute(sql) + # Get column names + rowid = {} + counter = 0 + for i in cursor.description: + rowid[i[0]] = counter + counter = counter + 1 + # print rowid + # Fetch all the rows in a list of lists. + results = cursor.fetchall() + for row in results: + result_dic[row["UniqueCruiseID"] + "_" + row["ConsecutiveCastNo"]] = { + keys: row[keys] for val, keys in enumerate(row.keys()) + } + return result_dic + except: + print("Error: unable to fecth data") + + +def read_mooring(db, cursor, table, MooringID): + sql = ("SELECT * from `{0}` WHERE `MooringID`= '{1}' ").format(table, MooringID) + + result_dic = {} + try: + # Execute the SQL command + cursor.execute(sql) + # Get column names + rowid = {} + counter = 0 + for i in cursor.description: + rowid[i[0]] = counter + counter = counter + 1 + # print rowid + # Fetch all the rows in a list of lists. + results = cursor.fetchall() + for row in results: + result_dic[row["MooringID"]] = { + keys: row[keys] for val, keys in enumerate(row.keys()) + } + return result_dic + except: + print("Error: unable to fecth data") + + +"""------------------------------------------------------------------------------------""" +parser = argparse.ArgumentParser( + description="Find Closest CTD casts to Mooring Deployment" +) + +parser.add_argument( + "DistanceThreshold", + metavar="DistanceThreshold", + type=float, + help="Distance From Mooring in Kilometers", +) +parser.add_argument( + "YearRange", + metavar="YearRange", + type=int, + nargs=2, + help="Range of years to look (eg 2012 2014)", +) +parser.add_argument( + "-db_moor", + "--db_moorings", + type=str, + help="path to db .pyini file for mooring records", + default='_secret/db_config_mooring.yaml' +) +parser.add_argument( + "-db_ctd", + "--db_ctd", + type=str, + help="path to db .pyini file for ctd records", + default='_secret/db_config_cruises.yaml' +) +parser.add_argument( + "-MooringID", metavar="--MooringID", type=str, help="MooringID 13BSM-2A" +) +parser.add_argument( + "-latlon", + "--latlon", + nargs="+", + type=float, + help="use manual lat/lon (decimaldegrees +N,+W)", +) + +args = parser.parse_args() + +host = "akutan" + +if not args.latlon and not args.MooringID: + print("Choose either a mooring location or a lat/lon pairing") + sys.exit() + +if args.latlon: # manual input of lat/lon + location = [args.latlon[0], args.latlon[1]] + +if args.MooringID: + # get db meta information for mooring + ### connect to DB + (db, cursor) = connect_to_DB(db_config_file=args.db_moorings) + table = "mooringdeploymentlogs" + Mooring_Meta = read_mooring(db, cursor, table, args.MooringID) + close_DB(db) + + # location = [71 + 13.413/60., 164 + 14.98/60.] + location = [ + float(Mooring_Meta[args.MooringID]["Latitude"].split()[0]) + + float(Mooring_Meta[args.MooringID]["Latitude"].split()[1]) / 60.0, + float(Mooring_Meta[args.MooringID]["Longitude"].split()[0]) + + float(Mooring_Meta[args.MooringID]["Longitude"].split()[1]) / 60.0, + ] + + +threshold = args.DistanceThreshold # km + +# get db meta information for mooring +### connect to DB +(db, cursor) = connect_to_DB(db_config_file=args.db_ctd) +table = "cruisecastlogs" +cruise_data = read_data(db, cursor, table, args.YearRange) +db.close() + +for index in sorted(cruise_data.keys()): + + destination = [ + cruise_data[index]["LatitudeDeg"] + cruise_data[index]["LatitudeMin"] / 60.0, + cruise_data[index]["LongitudeDeg"] + cruise_data[index]["LongitudeMin"] / 60.0, + ] + Distance2Station = sphered.distance(location, destination) + + if Distance2Station <= threshold: + print( + "Cast {0} on Cruise {1} is {2:3.2f} km away - {3}-{4}-{5} and {6}m deep".format( + cruise_data[index]["ConsecutiveCastNo"], + cruise_data[index]["UniqueCruiseID"], + Distance2Station, + cruise_data[index]["GMTYear"], + cruise_data[index]["GMTMonth"], + cruise_data[index]["GMTDay"], + cruise_data[index]["MaxDepth"], + ) + )