Patchgrid Improvement

Jenkinsfile

@@ -1,7 +1,14 @@
 
 pipeline {
     agent {
-        docker { image 'python:3.13' }
+        docker {
+            image 'python:3.13'
+            // The -u 0 flags means run commands inside the container
+            // as the user with uid = 0. This user is, by default, the
+            // root user. So it is effectively saying run the commands
+            // as root.
+            args '-u 0:0'
+        }
     }
     stages {
         stage('Installing OS Dependencies') {

qcore/grid.py

@@ -16,7 +16,6 @@
 from typing import Optional
 
 import numpy as np
-import scipy as sp
 
 from qcore import coordinates
 
@@ -204,11 +203,12 @@ def coordinate_patchgrid(
     origin: np.ndarray,
     x_upper: np.ndarray,
     y_bottom: np.ndarray,
-    resolution: float,
+    resolution: Optional[float] = None,
     nx: Optional[int] = None,
     ny: Optional[int] = None,
 ) -> np.ndarray:
-    """Creates a grid of patches in a bounded plane region.
+    """
+    Creates a grid of patches in a bounded plane region.
 
     Given the bounds of a rectangular planar region, create a grid of
     (lat, lon, depth) coordinates spaced at close to resolution metres apart
@@ -233,21 +233,51 @@ def coordinate_patchgrid(
     Returns
     -------
     np.ndarray
-        The patch grid of the rectangular planar region. Has shape (ny, nx), where
+        The patch grid of the rectangular planar region. Has shape (ny, nx, 3), where
         ny is the number of points in the origin->y_bottom direction and nx the number of
         points in the origin->x_upper direction.
+
+    Note
+    ----
+    The patch grid may have different sizes than given in as resolution if the resolution does not divide the lengths of the sides of the plane evenly.
+
+    Example
+    -------
+    >>> origin = np.array([-43.5321, 172.6362, 0.0])  # Christchurch, NZ
+    >>> x_upper = np.array([-43.5311, 172.6462, 0.0]) # ~800m to the east
+    >>> y_bottom = np.array([-43.5421, 172.6362, 0.0]) # ~1.2km to the south
+    >>> resolution = 100  # 100 meters
+    >>> grid = coordinate_patchgrid(origin, x_upper, y_bottom, resolution)
+    >>> grid.shape
+    (12, 8, 3)
     """
-    meshgrid = coordinate_meshgrid(origin, x_upper, y_bottom, resolution, nx=nx + 1 if nx is not None else nx, ny=ny + 1 if ny is not None else ny)
-    ny, nx = meshgrid.shape[:2]
-    meshgrid = coordinates.wgs_depth_to_nztm(meshgrid.reshape((-1, 3))).reshape(
-        (ny, nx, 3)
+    origin, x_upper, y_bottom = map(
+        coordinates.wgs_depth_to_nztm, (origin, x_upper, y_bottom)
     )
-    kernel = np.full((2, 2), 1 / 4)
-    patch_grid = np.zeros((ny - 1, nx - 1, 3))
-    for i in range(3):
-        patch_grid[:, :, i] = sp.signal.convolve2d(
-            meshgrid[:, :, i].reshape((ny, nx)), kernel, mode="valid"
-        )
+
+    v_x = x_upper - origin
+    v_y = y_bottom - origin
+
+    len_x = np.linalg.norm(v_x)
+    len_y = np.linalg.norm(v_y)
+
+    nx = nx or max(1, round(len_x / resolution))
+    ny = ny or max(1, round(len_y / resolution))
+
+    alpha, beta = np.meshgrid(
+        # The 1 / (2 * nx) term is to ensure that the patches are centred on the grid points.
+        # 1 / nx is the length of a patch, so the centre of the patch is 1 / (2 * nx) from the edge.
+        # the last patch is 1 / (2 * nx) from the edge, so the last grid point at 1 - 1 / (2 * nx).
+        # Similarly for ny.
+        np.linspace(1 / (2 * nx), 1 - 1 / (2 * nx), nx),
+        np.linspace(1 / (2 * ny), 1 - 1 / (2 * ny), ny),
+        indexing="ij",
+    )
+
+    patch_grid = np.transpose(
+        origin + alpha[..., np.newaxis] * v_x + beta[..., np.newaxis] * v_y, (1, 0, 2)
+    )
+
     return coordinates.nztm_to_wgs_depth(patch_grid.reshape((-1, 3))).reshape(
         patch_grid.shape
     )

qcore/test/test_grid/test_grid.py

@@ -1,7 +1,9 @@
-from qcore import grid, coordinates
-import pytest
 import itertools
+
 import numpy as np
+import pytest
+
+from qcore import coordinates, grid
 
 
 def test_grid_corners():
@@ -15,27 +17,43 @@ def test_grid_corners():
     projected_width = 5  # km
     dip = np.arctan((dbottom - dtop) / projected_width)
 
-    corners = grid.grid_corners(centroid, strike, dip_dir, dtop, dbottom, length, projected_width)
+    corners = grid.grid_corners(
+        centroid, strike, dip_dir, dtop, dbottom, length, projected_width
+    )
     assert corners.shape == (4, 3)  # Should return 4 corners with (lat, lon, depth)
-    assert corners[0, 2] == pytest.approx(dtop * 1000, abs=1e-2)  # Depth in meters for top corners
-    assert corners[2, 2] == pytest.approx(dbottom * 1000, abs=1e-2)  # Depth in meters for bottom corners
-    assert coordinates.distance_between_wgs_depth_coordinates(corners[1], corners[0]) == pytest.approx(length * 1000, abs=1e-2)
-    assert coordinates.distance_between_wgs_depth_coordinates(corners[2], corners[0]) == pytest.approx(projected_width / np.cos(dip) * 1000, abs=1e-2)
+    assert corners[0, 2] == pytest.approx(
+        dtop * 1000, abs=1e-2
+    )  # Depth in meters for top corners
+    assert corners[2, 2] == pytest.approx(
+        dbottom * 1000, abs=1e-2
+    )  # Depth in meters for bottom corners
+    assert coordinates.distance_between_wgs_depth_coordinates(
+        corners[1], corners[0]
+    ) == pytest.approx(length * 1000, abs=1e-2)
+    assert coordinates.distance_between_wgs_depth_coordinates(
+        corners[2], corners[0]
+    ) == pytest.approx(projected_width / np.cos(dip) * 1000, abs=1e-2)
     assert np.allclose(np.mean(corners, axis=0)[:2], centroid)
 
 
 def test_coordinate_meshgrid():
     """Test the coordinate_meshgrid function to verify grid creation"""
     origin = np.array([-43.5, 172.5, 5000])  # (lat, lon, depth in m)
-    x_upper = np.array([ -43.45498004,  172.50037108, 5000.        ])
-    y_bottom = np.array([ -43.50025372,  172.56184531, 5000.        ])
+    x_upper = np.array([-43.45498004, 172.50037108, 5000.0])
+    y_bottom = np.array([-43.50025372, 172.56184531, 5000.0])
     resolution = 1000  # 1 km resolution
 
     meshgrid = grid.coordinate_meshgrid(origin, x_upper, y_bottom, resolution)
     assert meshgrid.shape[2] == 3  # Should have shape (ny, nx, 3) for (lat, lon, depth)
-    for i, j in itertools.product(range(meshgrid.shape[0] - 1), range(meshgrid.shape[1]- 1) ):
-        assert coordinates.distance_between_wgs_depth_coordinates(meshgrid[i + 1, j], meshgrid[i, j]) == pytest.approx(resolution, abs=1e-2)
-        assert coordinates.distance_between_wgs_depth_coordinates(meshgrid[i, j + 1], meshgrid[i, j]) == pytest.approx(resolution, abs=1e-2)
+    for i, j in itertools.product(
+        range(meshgrid.shape[0] - 1), range(meshgrid.shape[1] - 1)
+    ):
+        assert coordinates.distance_between_wgs_depth_coordinates(
+            meshgrid[i + 1, j], meshgrid[i, j]
+        ) == pytest.approx(resolution, abs=1e-2)
+        assert coordinates.distance_between_wgs_depth_coordinates(
+            meshgrid[i, j + 1], meshgrid[i, j]
+        ) == pytest.approx(resolution, abs=1e-2)
 
     ny, nx = meshgrid.shape[:2]
     assert nx > 1 and ny > 1  # There should be multiple grid points in both directions
@@ -44,29 +62,75 @@ def test_coordinate_meshgrid():
     assert np.all(meshgrid[:, :, 2] == pytest.approx(origin[2], rel=1e-2))
 
 
+def test_coordinate_meshgrid_x_boundary():
+    """Test the coordinate_meshgrid function when you only have one patch in x and y direction"""
+    origin = np.array([-43.5, 172.5, 5000])  # (lat, lon, depth in m)
+    x_upper = np.array([-43.45498004, 172.50037108, 5000.0])
+    y_bottom = np.array([-43.50025372, 172.56184531, 5000.0])
+    resolution = coordinates.distance_between_wgs_depth_coordinates(
+        x_upper, origin
+    )  # 1 km resolution
+
+    meshgrid = grid.coordinate_patchgrid(origin, x_upper, y_bottom, resolution)
+    assert meshgrid.shape[2] == 3  # Should have shape (ny, nx, 3) for (lat, lon, depth)
+    ny, nx = meshgrid.shape[:2]
+    assert meshgrid.shape == (1, 1, 3)
+    # This one grid point should be the centre of the trapezium
+    origin_nztm = coordinates.wgs_depth_to_nztm(origin)
+    x_upper_nztm = coordinates.wgs_depth_to_nztm(x_upper)
+    y_bottom_nztm = coordinates.wgs_depth_to_nztm(y_bottom)
+    centre = x_upper_nztm + y_bottom_nztm - origin_nztm
+
+    assert meshgrid[0, 0] == pytest.approx(
+        coordinates.nztm_to_wgs_depth(centre), rel=1e-2
+    )
+
+
+def test_coordinate_patchgrid_different_nx_ny():
+    """Test the coordinate_meshgrid function to verify nx/ny parameters"""
+    origin = np.array([-43.5, 172.5, 5000])  # (lat, lon, depth in m)
+    x_upper = np.array([-43.45498004, 172.50037108, 5000.0])
+    y_bottom = np.array([-43.50025372, 172.56184531, 5000.0])
+    ny = 3
+    nx = 5
+    meshgrid = grid.coordinate_patchgrid(origin, x_upper, y_bottom, nx=nx, ny=ny)
+    assert meshgrid.shape == (ny, nx, 3)
+
+
 def test_coordinate_patchgrid():
     """Test the coordinate_meshgrid function to verify grid creation"""
     origin = np.array([-43.5, 172.5, 5000])  # (lat, lon, depth in m)
-    x_upper = np.array([ -43.45498004,  172.50037108, 5000.        ])
-    y_bottom = np.array([ -43.50025372,  172.56184531, 5000.        ])
+    x_upper = np.array([-43.45498004, 172.50037108, 5000.0])
+    y_bottom = np.array([-43.50025372, 172.56184531, 5000.0])
     resolution = 1000  # 1 km resolution
 
     meshgrid = grid.coordinate_patchgrid(origin, x_upper, y_bottom, resolution)
     assert meshgrid.shape[2] == 3  # Should have shape (ny, nx, 3) for (lat, lon, depth)
-    for i, j in itertools.product(range(meshgrid.shape[0] - 1), range(meshgrid.shape[1]- 1) ):
-        assert coordinates.distance_between_wgs_depth_coordinates(meshgrid[i + 1, j], meshgrid[i, j]) == pytest.approx(resolution, abs=1e-2)
-        assert coordinates.distance_between_wgs_depth_coordinates(meshgrid[i, j + 1], meshgrid[i, j]) == pytest.approx(resolution, abs=1e-2)
+    for i, j in itertools.product(
+        range(meshgrid.shape[0] - 1), range(meshgrid.shape[1] - 1)
+    ):
+        assert coordinates.distance_between_wgs_depth_coordinates(
+            meshgrid[i + 1, j], meshgrid[i, j]
+        ) == pytest.approx(resolution, abs=1e-2)
+        assert coordinates.distance_between_wgs_depth_coordinates(
+            meshgrid[i, j + 1], meshgrid[i, j]
+        ) == pytest.approx(resolution, abs=1e-2)
 
     ny, nx = meshgrid.shape[:2]
+    width = coordinates.distance_between_wgs_depth_coordinates(y_bottom, origin)
+    length = coordinates.distance_between_wgs_depth_coordinates(x_upper, origin)
+    assert ny == round(length / resolution)
+    assert nx == round(width / resolution)
+
     assert nx > 1 and ny > 1  # There should be multiple grid points in both directions
 
     # Check if the depths are consistent
     assert np.all(meshgrid[:, :, 2] == pytest.approx(origin[2], rel=1e-2))
 
 
 @pytest.mark.parametrize(
-    'length,resolution,expected',
-    [(10, 5, 3), (10, 3, 5),  (10000, 1000, 11), (1500, 500, 4)]
+    "length,resolution,expected",
+    [(10, 5, 3), (10, 3, 5), (10000, 1000, 11), (1500, 500, 4)],
 )
 def test_gridpoint_count_in_length(length: float, resolution: float, expected: int):
     """Test gridpoint_count_in_length function"""