From 159b6242622fd05412b908ccad89c9d11d97e173 Mon Sep 17 00:00:00 2001 From: Andy Beer Date: Tue, 26 Dec 2023 13:37:35 -0300 Subject: [PATCH 1/2] Crs handling conversion between different axis orderings --- resqpy/crs.py | 60 ++++++++++++- tests/unit_tests/test_crs.py | 163 +++++++++++++++++++++++++++++++++-- 2 files changed, 214 insertions(+), 9 deletions(-) diff --git a/resqpy/crs.py b/resqpy/crs.py index 0ed45ccb..d8366e41 100644 --- a/resqpy/crs.py +++ b/resqpy/crs.py @@ -22,6 +22,16 @@ PointType = Union[Tuple[float, float, float], List[float], np.ndarray] +axis_reordering = { # dictionary mapping from axis ordering to 2x2 transform matrix to convert from key to "easting northing" + "easting northing": np.array([(1.0, 0.0), (0.0, 1.0)], dtype = float), + "northing easting": np.array([(0.0, 1.0), (1.0, 0.0)], dtype = float), + "westing southing": np.array([(-1.0, 0.0), (0.0, -1.0)], dtype = float), + "southing westing": np.array([(0.0, -1.0), (-1.0, 0.0)], dtype = float), + "northing westing": np.array([(0.0, 1.0), (-1.0, 0.0)], dtype = float), + "westing northing": np.array([(-1.0, 0.0), (0.0, 1.0)], dtype = float), +} + + class Crs(BaseResqpy): """Coordinate reference system object.""" @@ -180,10 +190,17 @@ def is_right_handed_xyz(self) -> bool: return self.is_right_handed_xy() is bool(self.z_inc_down) - def global_to_local(self, xyz: PointType, global_z_inc_down: bool = True) -> Tuple[float, float, float]: + def global_to_local(self, + xyz: PointType, + global_z_inc_down: bool = True, + global_axis_order = 'easting northing') -> Tuple[float, float, float]: """Convert a single xyz point from the parent coordinate reference system to this one.""" x, y, z = xyz + # note: x & y offsets assumed to be in local axis ordering + if self.axis_order != global_axis_order: + assert global_axis_order in self.valid_axis_orders, f'invalid axis order: {global_axis_order}' + x, y = switch_axes(global_axis_order, self.axis_order, x, y) if self.x_offset != 0.0: x -= self.x_offset if self.y_offset != 0.0: @@ -196,9 +213,16 @@ def global_to_local(self, xyz: PointType, global_z_inc_down: bool = True) -> Tup (x, y, z) = vec.rotate_vector(self.rotation_matrix, np.array((x, y, z))) return (x, y, z) - def global_to_local_array(self, xyz: np.ndarray, global_z_inc_down: bool = True): + def global_to_local_array(self, + xyz: np.ndarray, + global_z_inc_down: bool = True, + global_axis_order = 'easting northing'): """Convert in situ a numpy array of xyz points from the parent coordinate reference system to this one.""" + # note: x & y offsets assumed to be in local axis ordering + if self.axis_order != global_axis_order: + assert global_axis_order in self.valid_axis_orders, f'invalid axis order: {global_axis_order}' + switch_axes_array(global_axis_order, self.axis_order, xyz) if self.x_offset != 0.0: xyz[..., 0] -= self.x_offset if self.y_offset != 0.0: @@ -212,7 +236,10 @@ def global_to_local_array(self, xyz: np.ndarray, global_z_inc_down: bool = True) a = vec.rotate_array(self.rotation_matrix, xyz) xyz[:] = a - def local_to_global(self, xyz: PointType, global_z_inc_down: bool = True) -> Tuple[float, float, float]: + def local_to_global(self, + xyz: PointType, + global_z_inc_down: bool = True, + global_axis_order = 'easting northing') -> Tuple[float, float, float]: """Convert a single xyz point from this coordinate reference system to the parent one.""" if self.rotated: @@ -227,9 +254,16 @@ def local_to_global(self, xyz: PointType, global_z_inc_down: bool = True) -> Tup z += self.z_offset if global_z_inc_down != self.z_inc_down: z = -z + # note: x & y offsets assumed to be in local axis ordering + if self.axis_order != global_axis_order: + assert global_axis_order in self.valid_axis_orders, f'invalid axis order: {global_axis_order}' + x, y = switch_axes(self.axis_order, global_axis_order, x, y) return (x, y, z) - def local_to_global_array(self, xyz: np.ndarray, global_z_inc_down: bool = True): + def local_to_global_array(self, + xyz: np.ndarray, + global_z_inc_down: bool = True, + global_axis_order = 'easting northing'): """Convert in situ a numpy array of xyz points from this coordinate reference system to the parent one.""" if self.rotated: @@ -244,6 +278,10 @@ def local_to_global_array(self, xyz: np.ndarray, global_z_inc_down: bool = True) if global_z_inc_down != self.z_inc_down: z = np.negative(xyz[..., 2]) xyz[..., 2] = z + # note: x & y offsets assumed to be in local axis ordering + if self.axis_order != global_axis_order: + assert global_axis_order in self.valid_axis_orders, f'invalid axis order: {global_axis_order}' + switch_axes_array(self.axis_order, global_axis_order, xyz) def has_same_epsg_code(self, other_crs: 'Crs') -> bool: """Returns True if either of the crs'es has a null EPSG code, or if they are the same.""" @@ -497,6 +535,20 @@ def create_xml(self, return crs +def switch_axes(from_order, to_order, x, y): + """Return x, y pair converting from one axis order to another.""" + matrix = np.matmul(axis_reordering[to_order].T, axis_reordering[from_order]) + return tuple(np.matmul(np.array((x, y), dtype = float), matrix)) + + +def switch_axes_array(from_order, to_order, xyz_array): + """Modify x & y values of array in situ, converting from one axis order to another.""" + matrix = np.matmul(axis_reordering[to_order].T, axis_reordering[from_order]) + assert xyz_array.shape[-1] in [2, 3] + xy = np.matmul(xyz_array[..., :2], matrix) + xyz_array[..., :2] = xy + + def _as_xyz_tuple(xyz): """Coerce into 3-tuple of floats.""" diff --git a/tests/unit_tests/test_crs.py b/tests/unit_tests/test_crs.py index b80de1f2..fb65df78 100644 --- a/tests/unit_tests/test_crs.py +++ b/tests/unit_tests/test_crs.py @@ -11,18 +11,150 @@ import resqpy.olio.uuid as bu +def test_switch_axes(): + xy = rqc.switch_axes("easting northing", "easting northing", 2.0, 3.0) + assert_array_almost_equal(xy, (2.0, 3.0)) + xy = rqc.switch_axes("northing easting", "northing easting", 2.0, 3.0) + assert_array_almost_equal(xy, (2.0, 3.0)) + xy = rqc.switch_axes("easting northing", "northing easting", 2.0, 3.0) + assert_array_almost_equal(xy, (3.0, 2.0)) + xy = rqc.switch_axes("northing easting", "easting northing", -2.0, 3.0) + assert_array_almost_equal(xy, (3.0, -2.0)) + xy = rqc.switch_axes("easting northing", "westing southing", 2.0, -3.0) + assert_array_almost_equal(xy, (-2.0, 3.0)) + xy = rqc.switch_axes("westing southing", "easting northing", 2.0, 3.0) + assert_array_almost_equal(xy, (-2.0, -3.0)) + xy = rqc.switch_axes("northing easting", "westing southing", -2.0, -3.0) + assert_array_almost_equal(xy, (3.0, 2.0)) + xy = rqc.switch_axes("westing southing", "northing easting", 2.0, 3.0) + assert_array_almost_equal(xy, (-3.0, -2.0)) + xy = rqc.switch_axes("westing southing", "northing easting", 0.0, 0.0) + assert_array_almost_equal(xy, (0.0, 0.0)) + xy = rqc.switch_axes("westing southing", "westing southing", 2.0, 3.0) + assert_array_almost_equal(xy, (2.0, 3.0)) + xy = rqc.switch_axes("easting northing", "southing westing", 2.0, 3.0) + assert_array_almost_equal(xy, (-3.0, -2.0)) + xy = rqc.switch_axes("southing westing", "northing easting", 2.0, 3.0) + assert_array_almost_equal(xy, (-2.0, -3.0)) + xy = rqc.switch_axes("easting northing", "northing westing", 2.0, 3.0) + assert_array_almost_equal(xy, (3.0, -2.0)) + xy = rqc.switch_axes("northing westing", "easting northing", 2.0, 3.0) + assert_array_almost_equal(xy, (-3.0, 2.0)) + xy = rqc.switch_axes("northing westing", "westing southing", 2.0, 3.0) + assert_array_almost_equal(xy, (3.0, -2.0)) + xy = rqc.switch_axes("westing southing", "northing westing", 2.0, 3.0) + assert_array_almost_equal(xy, (-3.0, 2.0)) + xy = rqc.switch_axes("westing northing", "westing southing", 2.0, 3.0) + assert_array_almost_equal(xy, (2.0, -3.0)) + xy = rqc.switch_axes("westing southing", "westing northing", 2.0, 3.0) + assert_array_almost_equal(xy, (2.0, -3.0)) + + +def test_switch_axes_array(): + a = np.array([(2.0, 3.0, 5.0), (7.0, 11.0, 13.0), (17.0, 19.0, 23.0), (29.0, 31.0, 37.0)], dtype = float) + xyz = a.copy() + rqc.switch_axes_array("easting northing", "easting northing", xyz) + assert_array_almost_equal(xyz, a) + xyz = a.copy() + rqc.switch_axes_array("northing easting", "northing easting", xyz) + assert_array_almost_equal(xyz, a) + xyz = a.copy() + rqc.switch_axes_array("easting northing", "northing easting", xyz) + assert_array_almost_equal(xyz[:, 0], a[:, 1]) + assert_array_almost_equal(xyz[:, 1], a[:, 0]) + assert_array_almost_equal(xyz[:, 2], a[:, 2]) + xyz = a.copy() + rqc.switch_axes_array("northing easting", "easting northing", xyz) + assert_array_almost_equal(xyz[:, 0], a[:, 1]) + assert_array_almost_equal(xyz[:, 1], a[:, 0]) + assert_array_almost_equal(xyz[:, 2], a[:, 2]) + xyz = a.copy() + rqc.switch_axes_array("easting northing", "westing southing", xyz) + assert_array_almost_equal(xyz[:, 0], -a[:, 0]) + assert_array_almost_equal(xyz[:, 1], -a[:, 1]) + assert_array_almost_equal(xyz[:, 2], a[:, 2]) + xyz = a.copy() + rqc.switch_axes_array("westing southing", "easting northing", xyz) + assert_array_almost_equal(xyz[:, 0], -a[:, 0]) + assert_array_almost_equal(xyz[:, 1], -a[:, 1]) + xyz = a.copy() + rqc.switch_axes_array("northing easting", "westing southing", xyz) + assert_array_almost_equal(xyz[:, 0], -a[:, 1]) + assert_array_almost_equal(xyz[:, 1], -a[:, 0]) + xy = a[:, :2].copy() # check working for xy data only + rqc.switch_axes_array("westing southing", "northing easting", xy) + assert_array_almost_equal(xy[:, 0], -a[:, 1]) + assert_array_almost_equal(xy[:, 1], -a[:, 0]) + xyz = a.copy() + rqc.switch_axes_array("westing southing", "northing easting", xyz) + assert_array_almost_equal(xyz[:, 0], -a[:, 1]) + assert_array_almost_equal(xyz[:, 1], -a[:, 0]) + xyz = a.copy() + rqc.switch_axes_array("westing southing", "westing southing", xyz) + assert_array_almost_equal(xyz, a) + xyz = a.copy() + rqc.switch_axes_array("easting northing", "southing westing", xyz) + assert_array_almost_equal(xyz[:, 0], -a[:, 1]) + assert_array_almost_equal(xyz[:, 1], -a[:, 0]) + xyz = a.copy() + rqc.switch_axes_array("southing westing", "northing easting", xyz) + assert_array_almost_equal(xyz[:, 0], -a[:, 0]) + assert_array_almost_equal(xyz[:, 1], -a[:, 1]) + xyz = a.copy() + rqc.switch_axes_array("easting northing", "northing westing", xyz) + assert_array_almost_equal(xyz[:, 0], a[:, 1]) + assert_array_almost_equal(xyz[:, 1], -a[:, 0]) + xyz = a.copy() + rqc.switch_axes_array("northing westing", "easting northing", xyz) + assert_array_almost_equal(xyz[:, 0], -a[:, 1]) + assert_array_almost_equal(xyz[:, 1], a[:, 0]) + xyz = a.copy() + rqc.switch_axes_array("northing westing", "westing southing", xyz) + assert_array_almost_equal(xyz[:, 0], a[:, 1]) + assert_array_almost_equal(xyz[:, 1], -a[:, 0]) + xyz = a.copy() + rqc.switch_axes_array("westing southing", "northing westing", xyz) + assert_array_almost_equal(xyz[:, 0], -a[:, 1]) + assert_array_almost_equal(xyz[:, 1], a[:, 0]) + xyz = a.copy() + rqc.switch_axes_array("westing northing", "westing southing", xyz) + assert_array_almost_equal(xyz[:, 0], a[:, 0]) + assert_array_almost_equal(xyz[:, 1], -a[:, 1]) + xyz = a.copy() + rqc.switch_axes_array("westing southing", "westing northing", xyz) + assert_array_almost_equal(xyz[:, 0], a[:, 0]) + assert_array_almost_equal(xyz[:, 1], -a[:, 1]) + # check also working for a single point + xyz = np.array((3.0, -5.0, 7.0), dtype = float) + rqc.switch_axes_array("northing westing", "westing southing", xyz) + assert_array_almost_equal(xyz, (-5.0, -3.0, 7.0)) + + def test_crs(tmp_path): # create some coordinate reference systems model = rq.new_model(os.path.join(tmp_path, 'crs_test.epc')) crs_default = rqc.Crs(model) assert crs_default.null_transform crs_m = rqc.Crs(model, xy_units = 'm', z_units = 'm') + crs_em_a = rqc.Crs(model, xy_units = 'm', z_units = 'm', extra_metadata = {'ab': 'a'}) + crs_em_b = rqc.Crs(model, xy_units = 'm', z_units = 'm', extra_metadata = {'ab': 'b'}) + crs_em_c = rqc.Crs(model, xy_units = 'm', z_units = 'm', extra_metadata = {'c': 'a'}) crs_ft = rqc.Crs(model, xy_units = 'ft', z_units = 'ft') crs_mixed = rqc.Crs(model, xy_units = 'm', z_units = 'ft') crs_offset = rqc.Crs(model, xy_units = 'm', z_units = 'm', x_offset = 100.0, y_offset = -100.0, z_offset = -50.0) assert not crs_offset.null_transform crs_elevation = rqc.Crs(model, z_inc_down = False) crs_rotate = rqc.Crs(model, rotation = maths.pi / 2.0, rotation_units = 'rad') + crs_north_rotate = rqc.Crs(model, + axis_order = 'northing easting', + rotation = maths.pi / 2.0, + rotation_units = 'rad') + crs_offset_rotate = rqc.Crs(model, + x_offset = 10.0, + y_offset = -10.0, + z_offset = -5.0, + rotation = maths.pi / 2.0, + rotation_units = 'rad') crs_south = rqc.Crs(model, axis_order = 'southing westing') crs_south_elevation = rqc.Crs(model, axis_order = 'southing westing', z_inc_down = False) crs_time_s = rqc.Crs(model, xy_units = 'm', time_units = 's') @@ -32,11 +164,16 @@ def test_crs(tmp_path): # check that distincitveness is recognised assert crs_default.is_equivalent(crs_m) + assert not crs_em_a == crs_m + assert not crs_em_b == crs_em_a + assert not crs_em_c == crs_em_a assert not crs_m.is_equivalent(crs_ft) assert not crs_mixed.is_equivalent(crs_m) assert not crs_m.is_equivalent(crs_offset) + assert not crs_m.is_equivalent(crs_offset_rotate) assert not crs_m.is_equivalent(crs_elevation) assert not crs_m.is_equivalent(crs_rotate) + assert not crs_m.is_equivalent(crs_north_rotate) assert not crs_m.is_equivalent(crs_south) assert not crs_time_s.is_equivalent(crs_time_ms) assert not crs_south_elevation.z_inc_down @@ -49,7 +186,9 @@ def test_crs(tmp_path): assert not crs_m.is_right_handed_xy() assert not crs_m.is_right_handed_xyz() assert not crs_elevation.is_right_handed_xy() + assert not crs_offset_rotate.is_right_handed_xy() assert crs_elevation.is_right_handed_xyz() + assert crs_north_rotate.is_right_handed_xy() assert crs_south.is_right_handed_xy() assert crs_south.is_right_handed_xyz() assert crs_south_elevation.is_right_handed_xy() @@ -58,7 +197,7 @@ def test_crs(tmp_path): # create some xml for crs in [ crs_default, crs_m, crs_ft, crs_mixed, crs_offset, crs_elevation, crs_rotate, crs_south, crs_time_s, - crs_time_ms + crs_time_ms, crs_offset_rotate, crs_em_a ]: crs.create_xml() model.store_epc() @@ -75,23 +214,37 @@ def test_crs(tmp_path): b = a.copy() crs_m.convert_array_from(crs_default, a) - assert np.max(np.abs(b - a)) < 1.0e-6 + assert_array_almost_equal(a, b) a[:] = b crs_m.convert_array_to(crs_m, a) assert np.all(a == b) crs_ft.convert_array_from(crs_m, a) - assert np.max(np.abs(b - a * ft_to_m)) < 1.0e-6 + assert_array_almost_equal(b, a * ft_to_m) crs_ft.convert_array_to(crs_m, a) - assert np.max(np.abs(b - a)) < 1.0e-6 + assert_array_almost_equal(a, b) a[:] = b crs_m.local_to_global_array(a) - assert np.max(np.abs(b - a)) < 1.0e-6 + assert_array_almost_equal(a, b) a[:] = b crs_offset.global_to_local_array(a) a[:, 0] += 100.0 a[:, 1] -= 100.0 a[:, 2] -= 50.0 assert_array_almost_equal(a, b) + a[:] = b + crs_mixed.convert_array_from(crs_m, a) + assert_array_almost_equal(b[..., :2], a[..., :2]) + assert_array_almost_equal(b[..., 2], a[..., 2] * ft_to_m) + a[:] = b + crs_south.convert_array_from(crs_m, a) + assert_array_almost_equal(a[..., 0], -b[..., 1]) + assert_array_almost_equal(a[..., 1], -b[..., 0]) + assert_array_almost_equal(a[..., 2], b[..., 2]) + a[:] = b + crs_south_elevation.convert_array_from(crs_m, a) + assert_array_almost_equal(a[..., 0], -b[..., 1]) + assert_array_almost_equal(a[..., 1], -b[..., 0]) + assert_array_almost_equal(a[..., 2], -b[..., 2]) # test single point conversion p = (456.78, 678.90, -1234.56) From cf3ea12d46cbcb2112eeb5febe7514b9192d7864 Mon Sep 17 00:00:00 2001 From: Andy Beer Date: Tue, 26 Dec 2023 15:29:43 -0300 Subject: [PATCH 2/2] unit tests --- tests/unit_tests/test_crs.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_crs.py b/tests/unit_tests/test_crs.py index fb65df78..d0a610cf 100644 --- a/tests/unit_tests/test_crs.py +++ b/tests/unit_tests/test_crs.py @@ -157,6 +157,7 @@ def test_crs(tmp_path): rotation_units = 'rad') crs_south = rqc.Crs(model, axis_order = 'southing westing') crs_south_elevation = rqc.Crs(model, axis_order = 'southing westing', z_inc_down = False) + crs_west = rqc.Crs(model, axis_order = 'westing southing') crs_time_s = rqc.Crs(model, xy_units = 'm', time_units = 's') crs_time_ms = rqc.Crs(model, xy_units = 'm', time_units = 'ms') for crs_time in [crs_time_s, crs_time_ms]: @@ -196,8 +197,8 @@ def test_crs(tmp_path): # create some xml for crs in [ - crs_default, crs_m, crs_ft, crs_mixed, crs_offset, crs_elevation, crs_rotate, crs_south, crs_time_s, - crs_time_ms, crs_offset_rotate, crs_em_a + crs_default, crs_m, crs_ft, crs_mixed, crs_offset, crs_elevation, crs_rotate, crs_south, crs_west, + crs_time_s, crs_time_ms, crs_offset_rotate, crs_em_a ]: crs.create_xml() model.store_epc() @@ -245,6 +246,18 @@ def test_crs(tmp_path): assert_array_almost_equal(a[..., 0], -b[..., 1]) assert_array_almost_equal(a[..., 1], -b[..., 0]) assert_array_almost_equal(a[..., 2], -b[..., 2]) + a[:] = b + crs_south.convert_array_to(crs_west, a) + assert_array_almost_equal(a[..., 0], b[..., 1]) + assert_array_almost_equal(a[..., 1], b[..., 0]) + assert_array_almost_equal(a[..., 2], b[..., 2]) + crs_south.convert_array_from(crs_west, a) + assert_array_almost_equal(a, b) + xyz = crs_south.convert_to(crs_west, a[0]) + assert len(xyz) == 3 + assert_array_almost_equal(xyz, (a[0, 1], a[0, 0], a[0, 2])) + xyz = crs_south.convert_from(crs_west, xyz) + assert_array_almost_equal(xyz, a[0]) # test single point conversion p = (456.78, 678.90, -1234.56)