Refactor crossings and lat/lon conversion functions

src/data/parse_mat.py

@@ -11,7 +11,7 @@ import pandas as pd
 from mat4py import loadmat
 from shapely.geometry import Point
 
-from data.parse_shp import convert_coord_systems
+from utils import convert_coord_systems
 from logs import setup_logging
 
 logger = setup_logging()

src/data/parse_shp.py

@@ -1,38 +1,20 @@
-from functools import partial
+"""
+This file can probably be removed since we are not reading from shp files any more
+"""
 
 import click
 import fiona
 import numpy as np
 import pandas as pd
-import pyproj
 from shapely.geometry import LineString, Point
 from shapely.geometry import shape
-from shapely.ops import transform
 
 from logs import setup_logging
+from utils import convert_coord_systems
 
 logger = setup_logging()
 
 
-def convert_coord_systems(g1, in_coord_system="EPSG:4326", out_coord_system="EPSG:28356"):
-    """
-    Converts coordinates from one coordinates system to another. Needed because shapefiles are usually defined in
-    lat/lon but should be converted to GDA to calculated distances.
-    https://gis.stackexchange.com/a/127432
-    :param in_coord_system: Default is lat/lon WGS84
-    :param out_coord_system: Default is GDA56 for NSW coastline
-    :return:
-    """
-    project = partial(
-        pyproj.transform,
-        pyproj.Proj(init=in_coord_system),  # source coordinate system
-        pyproj.Proj(init=out_coord_system),
-    )  # destination coordinate system
-
-    g2 = transform(project, g1)  # apply projection
-    return g2
-
-
 def shapes_from_shp(shp_file):
     """
     Parses a shape file and returns a list of shapely shapes, ids and properties

src/utils.py

@@ -1,15 +1,55 @@
-import os
-import yaml
+from functools import partial
 
+import numpy as np
+import pyproj
+from numpy import ma as ma
+from shapely.ops import transform
 
-def setup_logging(path="./src/logging.yaml", default_level=logging.INFO):
+
+def crossings(profile_x, profile_z, constant_z):
+    """
+    Finds the x coordinate of a z elevation for a beach profile. Much faster than using shapely to calculate
+    intersections since we are only interested in intersections of a constant value. Will return multiple
+    intersections if found. Used in calculating beach slope.
+    Adapted from https://stackoverflow.com/a/34745789
+    :param profile_x: List of x coordinates for the beach profile section
+    :param profile_z: List of z coordinates for the beach profile section
+    :param constant_z: Float of the elevation to find corresponding x coordinates
+    :return: List of x coordinates which correspond to the constant_z
+    """
+
+    # Remove nans to suppress warning messages
+    valid = ~ma.masked_invalid(profile_z).mask
+    profile_z = np.array(profile_z)[valid]
+    profile_x = np.array(profile_x)[valid]
+
+    # Normalize the 'signal' to zero.
+    # Use np.subtract rather than a list comprehension for performance reasons
+    z = np.subtract(profile_z, constant_z)
+
+    # Find all indices right before any crossing.
+    # TODO Sometimes this can give a runtime warning https://stackoverflow.com/a/36489085
+    indicies = np.where(z[:-1] * z[1:] < 0)[0]
+
+    # Use linear interpolation to find intersample crossings.
+    return [profile_x[i] - (profile_x[i] - profile_x[i + 1]) / (z[i] - z[i + 1]) * (z[i]) for i in indicies]
+
+
+# TODO Think of a better way to do this than having to manually specify the coordinate systems
+def convert_coord_systems(g1, in_coord_system="EPSG:4326", out_coord_system="EPSG:28356"):
     """
-    Setup logging configuration
+    Converts coordinates from one coordinates system to another. Needed because shapefiles are usually defined in
+    lat/lon but should be converted to GDA to calculated distances.
+    https://gis.stackexchange.com/a/127432
+    :param in_coord_system: Default is lat/lon WGS84
+    :param out_coord_system: Default is GDA56 for NSW coastline
+    :return:
     """
-    if os.path.exists(path):
-        with open(path, "rt") as f:
-            config = yaml.safe_load(f.read())
-        logging.config.dictConfig(config)
-    else:
-        logging.basicConfig(level=default_level)
-    return logging.getLogger(__name__)
+    project = partial(
+        pyproj.transform,
+        pyproj.Proj(init=in_coord_system),  # source coordinate system
+        pyproj.Proj(init=out_coord_system),
+    )  # destination coordinate system
+
+    g2 = transform(project, g1)  # apply projection
+    return g2