Skip to content

Commit

Permalink
work around for proj<9.1.0 2d-3d promotion
Browse files Browse the repository at this point in the history
  • Loading branch information
dugalh committed Oct 9, 2023
1 parent cb2d923 commit aaa9239
Show file tree
Hide file tree
Showing 3 changed files with 14 additions and 10 deletions.
4 changes: 3 additions & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,4 +28,6 @@
)
# TODO: the EPSG:<horiz>+<vert> format is not supported in rio 1.3.3, gdal 3.5.3, proj 9.1.0, but is supported
# in rio 1.3.6, gdal 3.6.2, proj 9.1.1. The exact version where support begins (proj=9.1.1?) should be set in
# setup.py. Then the multithreaded overview building bug should be fixed in gdal 3.8.
# setup.py (note that pyproj seems to support the EPSG:<>+<> format with PROJ<9.1.1 so something is exactly right
# with thre prev statement). Then the multithreaded overview building bug should be fixed in gdal 3.8. Also note
# that proj<9.1.0 auto promotes 2d CRSs to 3d.
9 changes: 4 additions & 5 deletions tests/test_io.py
Original file line number Diff line number Diff line change
Expand Up @@ -214,18 +214,17 @@ def test_rio_transform_vdatum_both(src_crs: str, dst_crs: str, request: pytest.F
('utm34n_crs', 'wgs84_egm96_crs'),
('utm34n_crs', 'wgs84_egm2008_crs'),
]) # yapf: disable
def _test_rio_transform_vdatum_one(src_crs: str, dst_crs: str, request: pytest.FixtureRequest):
def test_rio_transform_vdatum_one(src_crs: str, dst_crs: str, request: pytest.FixtureRequest):
"""
Test rasterio.warp.transform does not adjust the z coordinate with one of the source and destination CRS vertical
datums specified.
"""
src_crs: rio.CRS = rio.CRS.from_string(request.getfixturevalue(src_crs))
dst_crs: rio.CRS = rio.CRS.from_string(request.getfixturevalue(dst_crs))
src_xyz = [[10.], [10.], [100.]]
import os
os.environ.update(ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS='YES')
with rio.Env(ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS='YES'):
dst_xyz = transform(src_crs, dst_crs, *src_xyz)
dst_xyz = transform(src_crs, dst_crs, *src_xyz)
if rio.get_proj_version() >= (9, 1, 1):
# prior proj versions promote 2D->3D with ellipsoidal height
assert dst_xyz[2][0] == pytest.approx(src_xyz[2][0], abs=1e-6)


Expand Down
11 changes: 7 additions & 4 deletions tests/test_ortho.py
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ def test_reproject_dem_vdatum_both(
('float_utm34n_egm96_dem_file', 'utm34n_crs'),
],
) # yapf: disable # @formatter:on
def _test_reproject_dem_vdatum_one(
def test_reproject_dem_vdatum_one(
rgb_byte_src_file: Path, dem_file: str, pinhole_camera: Camera, crs: str, request: pytest.FixtureRequest
):
""" Test DEM reprojection does no altitude adjustment when one of DEM and ortho vertical datums are specified. """
Expand All @@ -234,7 +234,9 @@ def _test_reproject_dem_vdatum_one(
assert test_array.shape == ortho._dem_array.shape

mask = ~np.isnan(test_array) & ~np.isnan(ortho._dem_array)
assert test_array[mask] == pytest.approx(ortho._dem_array[mask], abs=1e-3)
if rio.get_proj_version() >= (9, 1, 0):
# prior proj versions promote 2D->3D with ellipsoidal height
assert test_array[mask] == pytest.approx(ortho._dem_array[mask], abs=1e-3)


@pytest.mark.parametrize('num_pts', [40, 100, 400, 1000, 4000])
Expand Down Expand Up @@ -316,8 +318,9 @@ def test_mask_dem(
assert dem_transform_mask == dem_transform
assert dem_mask.shape == ortho_mask.shape
assert dem_mask[ortho_mask].sum() / ortho_mask.sum() > 0.95
cc = np.corrcoef(dem_mask.flatten(), ortho_mask.flatten())
assert (np.all(dem_mask) and np.all(ortho_mask)) or (cc[0, 1] > 0.9)
if not np.all(dem_mask) and np.all(ortho_mask):
cc = np.corrcoef(dem_mask.flatten(), ortho_mask.flatten())
assert cc[0, 1] > 0.9

if False:
# debug plotting code
Expand Down

0 comments on commit aaa9239

Please sign in to comment.