Added transform_coordinates method to welleng.survey.SurveyParameters…

… for converting coordinates between projections. Added a couple of functions in utils for printing coordinates and better handling of conversions between decimal and dms coordinates.

tests/test_survey_parameters.py

@@ -1,24 +1,27 @@
 import unittest
 
 import welleng as we
+import numpy as np
 
 reference = {
     'x': 588319.02, 'y': 5770571.03, 'northing': 5770571.03,
     'easting': 588319.02, 'latitude': 52.077583926214494,
     'longitude': 4.288694821453205, 'convergence': 1.0166440347220762,
     'scale_factor': 0.9996957469340422, 'magnetic_field_intensity': 49381,
     'declination': 2.213, 'dip': -67.199, 'date': '2023-12-16',
-    'srs': 'EPSG:23031'
+    'srs': 'EPSG:23031',
+    'wgs84-utm31': [588225.162, 5770360.512]
 }
 
+calculator = we.survey.SurveyParameters(reference.get('srs'))
+
 
 class SurveyParamsTest(unittest.TestCase):
     def test_known_location(self):
-        calculator = we.survey.SurveyParameters(reference.get('srs'))
         survey_parameters = calculator.get_factors_from_x_y(
-            x=reference.get('x'), y=reference.get('y'), date=reference.get('date')
+            x=reference.get('x'), y=reference.get('y'),
+            date=reference.get('date')
         )
-
         for k, v in survey_parameters.items():
             try:
                 assert round(v, 3) == round(reference.get(k), 3)
@@ -27,6 +30,64 @@ def test_known_location(self):
 
         pass
 
+    def test_transform_projection_coordinates(self):
+        # Convert survey coordinates from UTM31_ED50 to UTM31_WGS84
+        coords = np.array((reference.get('easting'), reference.get('northing')))
+        result = calculator.transform_coordinates(coords, 'EPSG:32631')
+        assert np.allclose(
+            result,
+            np.array(reference.get('wgs84-utm31'))
+        )
+
+        # Try as a list
+        result = calculator.transform_coordinates(
+            coords.tolist(), 'EPSG:32631'
+        )
+        assert np.allclose(
+            result,
+            np.array(reference.get('wgs84-utm31'))
+        )
+
+        # Try as a tuple
+        result = calculator.transform_coordinates(
+            tuple(coords.tolist()), 'EPSG:32631'
+        )
+        assert np.allclose(
+            result,
+            np.array(reference.get('wgs84-utm31'))
+        )
+
+        result = calculator.transform_coordinates(
+            np.array([coords, coords]),
+            'EPSG:32631'
+        )
+        assert np.allclose(
+            result,
+            np.full_like(result, reference.get('wgs84-utm31'))
+        )
+
+        # Try as a list
+        result = calculator.transform_coordinates(
+            [coords.tolist(), coords.tolist()],
+            'EPSG:32631'
+        )
+        assert np.allclose(
+            result,
+            np.full_like(result, reference.get('wgs84-utm31'))
+        )
+
+        # Try as a tuple
+        result = calculator.transform_coordinates(
+            (tuple(coords.tolist()), tuple(coords.tolist())),
+            'EPSG:32631'
+        )
+        assert np.allclose(
+            result,
+            np.full_like(result, reference.get('wgs84-utm31'))
+        )
+
+        pass
+
 
 # def one_function_to_run_them_all():
 #     """

tests/test_utils.py

@@ -4,7 +4,9 @@
 from numpy.typing import NDArray
 
 from welleng.units import ureg
-from welleng.utils import annular_volume, decimal2dms, dms2decimal
+from welleng.utils import (
+    annular_volume, decimal2dms, dms2decimal, pprint_dms, dms_from_string
+)
 
 LAT, LON = (52, 4, 43.1868, 'N'), (4, 17, 19.6368, 'E')
 
@@ -40,6 +42,10 @@ def _generate_random_dms(n: int, ndigits: int = None) -> NDArray:
     ).reshape((-1, 2, 4))
 
 
+def are_tuples_identical(tuple1, tuple2):
+    return all(x == y for x, y in zip(tuple1, tuple2))
+
+
 class UtilsTest(unittest.TestCase):
     def test_annular_volume(self):
         av = annular_volume(
@@ -135,6 +141,29 @@ def test_dms2decimal2dms(self):
 
         assert np.all(np.equal(_dms, dms))
 
+    def test_pprint_dms(self):
+        result = pprint_dms(LAT, return_data=True)
+        data = dms_from_string(result)
+
+        assert are_tuples_identical(LAT, data)
+
+        result = pprint_dms(LAT, return_data=True, symbols=False)
+        data = dms_from_string(result)
+
+        assert are_tuples_identical(LAT, data)
+
+        result = pprint_dms(LAT[:3], return_data=True)
+        data = dms_from_string(result)
+
+        assert are_tuples_identical(LAT[:3], data)
+
+        result = pprint_dms(LAT[:3], return_data=True, symbols=False)
+        data = dms_from_string(result)
+
+        assert are_tuples_identical(LAT[:3], data)
+
+        pass
+
 
 # def one_function_to_run_them_all():
 #     """

welleng/survey.py

@@ -8,7 +8,8 @@
 except ImportError:
     MAG_CALC = False
 from datetime import datetime
-from pyproj import CRS, Proj
+from pyproj import CRS, Proj, Transformer
+from pyproj.enums import TransformDirection
 from scipy.optimize import minimize
 from scipy.spatial.transform import Rotation as R
 
@@ -31,7 +32,7 @@
 from .units import ureg
 
 from typing import List, Union
-from numpy.typing import NDArray
+from numpy.typing import NDArray, ArrayLike
 
 
 AZI_REF = ["true", "magnetic", "grid"]
@@ -48,6 +49,19 @@ class SurveyParameters(Proj):
     Requires ``pyproj`` and ``magnetic_field_calculator`` to be installed and
     access to the internet.
 
+    For reference, here's some EPSG codes:
+    {
+        'UTM31_ED50': 'EPSG:23031',
+        'UTM31_WGS84': 'EPSG:32631',
+        'RD': 'EPSG:28992',
+        'ED50-UTM31': 'EPSG:23031',
+        'ED50-NEDTM': 'EPSG:23095',  # assume same as ED50-UTM31
+        'ETRS89-UTM31': 'EPSG:25831',
+        'ED50-UTM32': 'EPSG:23032',
+        'ED50-GEOGR': 'EPSG:4230',
+        'WGS84-UTM31': 'EPSG:32631'
+    }
+
     References
     ----------
     For more info on transformations between maps, refer to the pyproj project
@@ -189,6 +203,53 @@ def get_factors_from_x_y(
 
         return data
 
+    def transform_coordinates(
+        self, coords: ArrayLike, to_projection: str, altitude: float = None,
+        *args, **kwargs
+    ) -> ArrayLike:
+        """Transforms coordinates from instance's projection to another
+        projection.
+
+        Parameters
+        ----------
+        coords: arraylike
+            A list of decimal coordinates to transform from the instance
+            projection to the specified projection system. Can be 2D or 3D in
+            (x, y, z) format, where x is East/West and y is North/South.
+        to_projection: str
+            The EPSG code of the desired coordinates.
+
+        Returns
+        -------
+        result: ArrayLike
+            An array of transformed coordinates in the desired projection.
+
+        Examples
+        --------
+        Convert the coordinates of Den Haag from ED50-UTM31 to WGS84-UTM31:
+
+        >>> from welleng.survey import SurveyParameters
+        >>> calculator = SurveyParameters('EPSG:23031')
+        >>> result = calculator.transform_coordinates(
+        ...     coords=[(588319.02, 5770571.03)], to_projection='EPSG:32631'
+        ... )
+        >>> print(result)
+        [[ 588225.93417027 5770360.56500115]]
+        """
+        transformer = Transformer.from_crs(
+            self.crs, CRS(to_projection)
+        )
+        _coords = np.array(coords)
+        result = list(transformer.itransform(
+            (
+                _coords.tolist() if len(_coords.shape) > 1
+                else _coords.reshape((1, -1)).tolist()
+            ),
+            direction=TransformDirection('FORWARD'), *args, **kwargs
+        ))
+
+        return np.array(result)
+
 
 class SurveyHeader:
     """

welleng/utils.py

@@ -1,3 +1,4 @@
+import re
 from typing import Annotated, Literal, Union
 
 import numpy as np
@@ -788,3 +789,57 @@ def dms2decimal(dms: tuple | NDArray, ndigits: int = None) -> NDArray:
         return result.reshape(-1)
     else:
         return result
+
+
+def pprint_dms(dms, symbols: bool = True, return_data: bool = False):
+    """Pretty prints a (decimal, minutes, seconds) tuple or list.
+
+    Parameters
+    ----------
+    dms: tuple | list
+        An x or y or northing or easting (degree, minute, second).
+    symbols: bool (default: True)
+        Whether to print symbols for (deg, min, sec).
+    return_data: bool (default: False)
+        If True then will return the string rather than print it.
+    """
+    if symbols:
+        try:
+            deg, min, sec = dms
+            text = f"{deg}\N{DEGREE SIGN}, {min}', {sec}\""
+        except ValueError:
+            deg, min, sec, _ = dms
+            text = f"{deg}\N{DEGREE SIGN}, {min}', {sec}\" {_}"
+
+    else:
+        try:
+            deg, min, sec = dms
+            text = f"{deg} deg, {min} min, {sec} sec"
+        except ValueError:
+            deg, min, sec, _ = dms
+            text = f"{deg} deg, {min} min, {sec} sec {_}"
+
+    if return_data:
+        return text
+    else:
+        print(text)
+
+
+def dms_from_string(text):
+    """Extracts the values from a string dms x or y or northing or easting.
+    """
+    pattern = re.compile(r'(\d+)\s*(?:°|deg)?,\s*(\d+)\s*(?:\'|min)?,\s*(\d+(?:\.\d+)?)\s*(sec)?\s*.*?(\S+)?$', re.IGNORECASE)
+    matches = pattern.findall(text)
+
+    if matches:
+        deg, min, sec_str = matches[0][:3]
+        sec = float(sec_str)
+        final_data = matches[0][-1] if matches[0][-1] else None
+
+        if final_data:
+            return (int(deg), int(min), sec, final_data)
+        else:
+            return (int(deg), int(min), sec)
+
+    else:
+        return