feat: add optional reference projection for CoordinateTransformer

seismostats/analysis/bvalue/base.py

@@ -44,7 +44,8 @@ def __call__(self,
         if not self.weights_supported and weights is not None:
             raise ValueError('Weights are not supported by this estimator.')
 
-        self.magnitudes = magnitudes
+        self.magnitudes = magnitudes.copy()
+
         self.weights = weights
 
         self._sanity_checks()

seismostats/analysis/bvalue/tests/test_base.py

@@ -1,5 +1,6 @@
 import warnings
 
+import numpy as np
 import pytest
 
 from seismostats.analysis.bvalue import ClassicBValueEstimator
@@ -19,3 +20,11 @@ def test_estimate_b_warnings():
         estimator = ClassicBValueEstimator(mc=-1, delta_m=0.1)
         estimator(mags)
         assert w[-1].category == UserWarning
+
+
+def test_by_reference():
+    mags = simulate_magnitudes_binned(n=100, b=1, mc=0, delta_m=0.1)
+    estimator = ClassicBValueEstimator(mc=1, delta_m=0.1)
+    estimator(mags)
+    estimator.magnitudes.sort()
+    assert not np.array_equal(mags, estimator.magnitudes)

seismostats/utils/coordinates.py

@@ -16,6 +16,12 @@ class CoordinateTransformer:
     "+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
     to represent a UTM coordinate system.
 
+    It is possible to pass a reference point, which will be used as the
+    origin of the local coordinate system. If no reference point is passed,
+    the origin will be set to (0, 0, 0). The projection of the reference point
+    is assumed to be the same as the local projection, unless specified
+    otherwise in the ``ref_proj`` kwarg.
+
     Notes:
         4326 as well as eg 2056 are 2D coordinate systems, so altitude
         is not taken into account and only calculated in reference to the ref.
@@ -28,23 +34,33 @@ def __init__(
             ref_easting: float = 0.0,
             ref_northing: float = 0.0,
             ref_altitude: float = 0.0,
-            external_proj: int | str = 4326):
+            external_proj: int | str = 4326,
+            ref_proj: int | str | None = None):
         """
         Constructor of CoordinateTransformer object.
 
         Args:
-            local_proj:     int (epsg) or string (proj) of local CRS.
-            ref_easting:    reference easting for local coordinates.
-            ref_northing:   reference northing for local coordinates.
-            ref_altitude:   reference altitude for local coordinates.
-            external_proj:  int or string of geographic coordinates.
+            local_proj:     projection of local CRS,
+                            eg. str('epsg:2056') or int(2056).
+            ref_easting:    reference easting
+            ref_northing:   reference northing
+            ref_altitude:   reference altitude
+            external_proj:  projection of external CRS,
+                            eg. str('epsg:4326') or int(4326).
+            ref_proj:       projection of reference coordinates,
+                            defaults to 'local_proj'.
         """
         self.ref_easting = ref_easting
         self.ref_northing = ref_northing
         self.ref_altitude = ref_altitude
         self.local_proj = local_proj
         self.external_proj = external_proj
 
+        if ref_proj is not None:
+            tr = Transformer.from_proj(ref_proj, local_proj, always_xy=True)
+            self.ref_easting, self.ref_northing = \
+                tr.transform(self.ref_easting, self.ref_northing)
+
         self.transformer_to_local = Transformer.from_proj(
             self.external_proj, self.local_proj, always_xy=True)
         self.transformer_to_external = Transformer.from_proj(

seismostats/utils/tests/test_coordinate_transformer.py

@@ -15,6 +15,9 @@
 REF_NORTHING = 1239374.8568319362
 REF_ALTITUDE = 0.0
 
+REF_LON = 8.54
+REF_LAT = 47.3
+
 # Coordinates for Bedretto lab
 BEDRETTO_LAB_ORIGIN_EASTING = 2679720.70
 BEDRETTO_LAB_ORIGIN_NORTHING = 1151600.13
@@ -94,3 +97,25 @@ def test_bedretto_transformation(self):
         self.assertAlmostEqual(lat, MOCK_INJ_LAT)
         self.assertAlmostEqual(lon, MOCK_INJ_LON)
         self.assertAlmostEqual(altitude, MOCK_INJ_ALTITUDE)
+
+    def test_reference_proj(self):
+        """
+        Test transformer with a reference point and a different
+        reference projection. The reference projection is set to
+        WGS84, so the reference point is transformed to the local
+        coordinate system.
+        """
+        transformer = CoordinateTransformer(
+            SWISS_PROJ, ref_easting=REF_LON, ref_northing=REF_LAT,
+            ref_proj=WGS84_PROJ)
+
+        easting, northing, altitude = transformer.\
+            to_local_coords(
+                LON, LAT)
+        self.assertAlmostEqual(easting, 0.0)
+        self.assertAlmostEqual(northing, 0.0)
+
+        lon, lat, _ = transformer.from_local_coords(
+            easting, northing)
+        self.assertAlmostEqual(lat, LAT)
+        self.assertAlmostEqual(lon, LON)