Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Crs handling conversion between different axis orderings #774

Merged
merged 2 commits into from
Dec 26, 2023
Merged

Crs handling conversion between different axis orderings #774

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
159b624
Crs handling conversion between different axis orderings
andy-beer Dec 26, 2023
cf3ea12
unit tests
andy-beer Dec 26, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 56 additions & 4 deletions resqpy/crs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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."""

Expand Down Expand Up @@ -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:
Expand All @@ -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:
Expand All @@ -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:
Expand All @@ -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:
Expand All @@ -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."""
Expand Down Expand Up @@ -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."""

Expand Down
178 changes: 172 additions & 6 deletions tests/unit_tests/test_crs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,32 +11,170 @@
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_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]:
assert crs_time.resqml_type == 'LocalTime3dCrs'

# 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
Expand All @@ -49,16 +187,18 @@ 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()
assert not crs_south_elevation.is_right_handed_xyz()

# 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_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()
Expand All @@ -75,23 +215,49 @@ 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])
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)
Expand Down