diff --git a/.gitignore b/.gitignore index 10b3931c..24590c58 100644 --- a/.gitignore +++ b/.gitignore @@ -146,6 +146,7 @@ xdem/_version.py # Example data downloaded/produced during tests examples/data/ +tests/test_data/ doc/source/basic_examples/ doc/source/advanced_examples/ diff --git a/Makefile b/Makefile index 21220f1d..8b9e66c9 100644 --- a/Makefile +++ b/Makefile @@ -11,17 +11,18 @@ ifndef VENV VENV = "venv" endif -# Python version requirement -PYTHON_VERSION_REQUIRED = 3.10 - +# Python global variables definition +PYTHON_VERSION_MIN = 3.10 +# Set PYTHON if not defined in command line +# Example: PYTHON="python3.10" make venv to use python 3.10 for the venv +# By default the default python3 of the system. ifndef PYTHON - # Try to find python version required - PYTHON = "python$(PYTHON_VERSION_REQUIRED)" + PYTHON = "python3" endif PYTHON_CMD=$(shell command -v $(PYTHON)) -PYTHON_VERSION_CUR=$(shell $(PYTHON_CMD) -c 'import sys; print("%d.%d" % sys.version_info[0:2])') -PYTHON_VERSION_OK=$(shell $(PYTHON_CMD) -c 'import sys; req_ver = tuple(map(int, "$(PYTHON_VERSION_REQUIRED)".split("."))); cur_ver = sys.version_info[0:2]; print(int(cur_ver == req_ver))') +PYTHON_VERSION_CUR=$(shell $(PYTHON_CMD) -c 'import sys; print("%d.%d"% sys.version_info[0:2])') +PYTHON_VERSION_OK=$(shell $(PYTHON_CMD) -c 'import sys; cur_ver = sys.version_info[0:2]; min_ver = tuple(map(int, "$(PYTHON_VERSION_MIN)".split("."))); print(int(cur_ver >= min_ver))') ############### Check python version supported ############ @@ -30,7 +31,7 @@ ifeq (, $(PYTHON_CMD)) endif ifeq ($(PYTHON_VERSION_OK), 0) - $(error "Requires Python version == $(PYTHON_VERSION_REQUIRED). Current version is $(PYTHON_VERSION_CUR)") + $(error "Requires Python version >= $(PYTHON_VERSION_MIN). Current version is $(PYTHON_VERSION_CUR)") endif ################ MAKE Targets ###################### @@ -45,19 +46,6 @@ venv: ## Create a virtual environment in 'venv' directory if it doesn't exist @touch ${VENV}/bin/activate @${VENV}/bin/python -m pip install --upgrade wheel setuptools pip -.PHONY: install-gdal -install-gdal: ## Install GDAL version matching the system's GDAL via pip - @if command -v gdalinfo >/dev/null 2>&1; then \ - GDAL_VERSION=$$(gdalinfo --version | awk '{print $$2}'); \ - echo "System GDAL version: $$GDAL_VERSION"; \ - ${VENV}/bin/pip install gdal==$$GDAL_VERSION; \ - else \ - echo "Warning: GDAL not found on the system. Proceeding without GDAL."; \ - echo "Try installing GDAL by running the following commands depending on your system:"; \ - echo "Debian/Ubuntu: sudo apt-get install -y gdal-bin libgdal-dev"; \ - echo "Red Hat/CentOS: sudo yum install -y gdal gdal-devel"; \ - echo "Then run 'make install-gdal' to proceed with GDAL installation."; \ - fi .PHONY: install install: venv ## Install xDEM for development (depends on venv) @@ -66,8 +54,6 @@ install: venv ## Install xDEM for development (depends on venv) @test -f .git/hooks/pre-commit || echo "Installing pre-commit hooks" @test -f .git/hooks/pre-commit || ${VENV}/bin/pre-commit install -t pre-commit @test -f .git/hooks/pre-push || ${VENV}/bin/pre-commit install -t pre-push - @echo "Attempting to install GDAL..." - @make install-gdal @echo "xdem installed in development mode in virtualenv ${VENV}" @echo "To use: source ${VENV}/bin/activate; xdem -h" diff --git a/dev-environment.yml b/dev-environment.yml index 73eddcd5..29e31bb2 100644 --- a/dev-environment.yml +++ b/dev-environment.yml @@ -9,7 +9,7 @@ dependencies: - matplotlib=3.* - pyproj>=3.4,<4 - rasterio>=1.3,<2 - - scipy=1.* + - scipy<1.15.0 - tqdm - scikit-image=0.* - scikit-gstat>=1.0.18,<1.1 @@ -28,13 +28,11 @@ dependencies: - scikit-learn # Test dependencies - - gdal # To test against GDAL - pytest - pytest-xdist - pyyaml - flake8 - pylint - - richdem # To test against richdem # Doc dependencies - sphinx diff --git a/environment.yml b/environment.yml index 3fe62ae4..153b68f9 100644 --- a/environment.yml +++ b/environment.yml @@ -9,7 +9,7 @@ dependencies: - matplotlib=3.* - pyproj>=3.4,<4 - rasterio>=1.3,<2 - - scipy=1.* + - scipy<1.15.0 - tqdm - scikit-image=0.* - scikit-gstat>=1.0.18,<1.1 diff --git a/requirements.txt b/requirements.txt index c293c1be..ca541404 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,7 +7,7 @@ numpy==1.* matplotlib==3.* pyproj>=3.4,<4 rasterio>=1.3,<2 -scipy==1.* +scipy<1.15.0 tqdm scikit-image==0.* scikit-gstat>=1.0.18,<1.1 diff --git a/setup.cfg b/setup.cfg index 8adfa7ae..45525c54 100644 --- a/setup.cfg +++ b/setup.cfg @@ -61,7 +61,6 @@ test = flake8 pylint scikit-learn - richdem doc = sphinx sphinx-book-theme diff --git a/tests/conftest.py b/tests/conftest.py index 675084e1..55ad79fb 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,142 +1,76 @@ -from typing import Callable, List, Union +import os +import tarfile +import tempfile +import urllib +from distutils.dir_util import copy_tree +from typing import Callable -import geoutils as gu -import numpy as np import pytest -import richdem as rd -from geoutils.raster import RasterType -from xdem._typing import NDArrayf +_TESTDATA_DIRECTORY = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "tests", "test_data")) +# Define a URL to the xdem-data repository's test data +_TESTDATA_REPO_URL = "https://github.com/vschaffn/xdem-data/tarball/2-richdem_gdal" +_COMMIT_HASH = "31a7159c982cec4b352f0de82bd4e0be61db3afe" -@pytest.fixture(scope="session") # type: ignore -def raster_to_rda() -> Callable[[RasterType], rd.rdarray]: - def _raster_to_rda(rst: RasterType) -> rd.rdarray: - """ - Convert geoutils.Raster to richDEM rdarray. - """ - arr = rst.data.filled(rst.nodata).squeeze() - rda = rd.rdarray(arr, no_data=rst.nodata) - rda.geotransform = rst.transform.to_gdal() - return rda - return _raster_to_rda +def download_test_data(overwrite: bool = False) -> None: + """ + Download the entire test_data directory from the xdem-data repository. + :param overwrite: If True, re-downloads the data even if it already exists. + """ + if not overwrite and os.path.exists(_TESTDATA_DIRECTORY) and os.listdir(_TESTDATA_DIRECTORY): + return # Test data already exists -@pytest.fixture(scope="session") # type: ignore -def get_terrainattr_richdem(raster_to_rda: Callable[[RasterType], rd.rdarray]) -> Callable[[RasterType, str], NDArrayf]: - def _get_terrainattr_richdem(rst: RasterType, attribute: str = "slope_radians") -> NDArrayf: - """ - Derive terrain attribute for DEM opened with geoutils.Raster using RichDEM. - """ - rda = raster_to_rda(rst) - terrattr = rd.TerrainAttribute(rda, attrib=attribute) - terrattr[terrattr == terrattr.no_data] = np.nan - return np.array(terrattr) + # Clear the directory if overwrite is True + if overwrite and os.path.exists(_TESTDATA_DIRECTORY): + for root, dirs, files in os.walk(_TESTDATA_DIRECTORY, topdown=False): + for name in files: + os.remove(os.path.join(root, name)) + for name in dirs: + os.rmdir(os.path.join(root, name)) + + # Create a temporary directory to download the tarball + temp_dir = tempfile.TemporaryDirectory() + tar_path = os.path.join(temp_dir.name, "test_data.tar.gz") + + # Construct the URL with the commit hash + url = f"{_TESTDATA_REPO_URL}#commit={_COMMIT_HASH}" + + response = urllib.request.urlopen(url) + if response.getcode() == 200: + with open(tar_path, "wb") as outfile: + outfile.write(response.read()) + else: + raise ValueError(f"Failed to download test data: {response.status_code}") - return _get_terrainattr_richdem + # Extract the tarball + with tarfile.open(tar_path) as tar: + tar.extractall(temp_dir.name) + + # Copy the test_data directory to the target directory + extracted_dir = os.path.join( + temp_dir.name, + [dirname for dirname in os.listdir(temp_dir.name) if os.path.isdir(os.path.join(temp_dir.name, dirname))][0], + "test_data", + ) + + copy_tree(extracted_dir, _TESTDATA_DIRECTORY) @pytest.fixture(scope="session") # type: ignore -def get_terrain_attribute_richdem( - get_terrainattr_richdem: Callable[[RasterType, str], NDArrayf] -) -> Callable[[RasterType, Union[str, list[str]], bool, float, float, float], Union[RasterType, list[RasterType]]]: - def _get_terrain_attribute_richdem( - dem: RasterType, - attribute: Union[str, List[str]], - degrees: bool = True, - hillshade_altitude: float = 45.0, - hillshade_azimuth: float = 315.0, - hillshade_z_factor: float = 1.0, - ) -> Union[RasterType, List[RasterType]]: - """ - Derive one or multiple terrain attributes from a DEM using RichDEM. - """ - if isinstance(attribute, str): - attribute = [attribute] - - if not isinstance(dem, gu.Raster): - raise ValueError("DEM must be a geoutils.Raster object.") - - terrain_attributes = {} - - # Check which products should be made to optimize the processing - make_aspect = any(attr in attribute for attr in ["aspect", "hillshade"]) - make_slope = any( - attr in attribute - for attr in [ - "slope", - "hillshade", - "planform_curvature", - "aspect", - "profile_curvature", - "maximum_curvature", - ] - ) - make_hillshade = "hillshade" in attribute - make_curvature = "curvature" in attribute - make_planform_curvature = "planform_curvature" in attribute or "maximum_curvature" in attribute - make_profile_curvature = "profile_curvature" in attribute or "maximum_curvature" in attribute - - if make_slope: - terrain_attributes["slope"] = get_terrainattr_richdem(dem, "slope_radians") - - if make_aspect: - # The aspect of RichDEM is returned in degrees, we convert to radians to match the others - terrain_attributes["aspect"] = np.deg2rad(get_terrainattr_richdem(dem, "aspect")) - # For flat slopes, RichDEM returns a 90° aspect by default, while GDAL return a 180° aspect - # We stay consistent with GDAL - slope_tmp = get_terrainattr_richdem(dem, "slope_radians") - terrain_attributes["aspect"][slope_tmp == 0] = np.pi - - if make_hillshade: - # If a different z-factor was given, slopemap with exaggerated gradients. - if hillshade_z_factor != 1.0: - slopemap = np.arctan(np.tan(terrain_attributes["slope"]) * hillshade_z_factor) - else: - slopemap = terrain_attributes["slope"] - - azimuth_rad = np.deg2rad(360 - hillshade_azimuth) - altitude_rad = np.deg2rad(hillshade_altitude) - - # The operation below yielded the closest hillshade to GDAL (multiplying by 255 did not work) - # As 0 is generally no data for this uint8, we add 1 and then 0.5 for the rounding to occur between - # 1 and 255 - terrain_attributes["hillshade"] = np.clip( - 1.5 - + 254 - * ( - np.sin(altitude_rad) * np.cos(slopemap) - + np.cos(altitude_rad) * np.sin(slopemap) * np.sin(azimuth_rad - terrain_attributes["aspect"]) - ), - 0, - 255, - ).astype("float32") - - if make_curvature: - terrain_attributes["curvature"] = get_terrainattr_richdem(dem, "curvature") - - if make_planform_curvature: - terrain_attributes["planform_curvature"] = get_terrainattr_richdem(dem, "planform_curvature") - - if make_profile_curvature: - terrain_attributes["profile_curvature"] = get_terrainattr_richdem(dem, "profile_curvature") - - # Convert the unit if wanted. - if degrees: - for attr in ["slope", "aspect"]: - if attr not in terrain_attributes: - continue - terrain_attributes[attr] = np.rad2deg(terrain_attributes[attr]) - - output_attributes = [terrain_attributes[key].reshape(dem.shape) for key in attribute] - - if isinstance(dem, gu.Raster): - output_attributes = [ - gu.Raster.from_array(attr, transform=dem.transform, crs=dem.crs, nodata=-99999) - for attr in output_attributes - ] - - return output_attributes if len(output_attributes) > 1 else output_attributes[0] - - return _get_terrain_attribute_richdem +def get_test_data_path() -> Callable[[str], str]: + def _get_test_data_path(filename: str, overwrite: bool = False) -> str: + """Get file from test_data""" + download_test_data(overwrite=overwrite) # Ensure the test data is downloaded + file_path = os.path.join(_TESTDATA_DIRECTORY, filename) + + if not os.path.exists(file_path): + if overwrite: + raise FileNotFoundError(f"The file {filename} was not found in the test_data directory.") + file_path = _get_test_data_path(filename, overwrite=True) + + return file_path + + return _get_test_data_path diff --git a/tests/test_coreg/test_affine.py b/tests/test_coreg/test_affine.py index a642038c..185398c1 100644 --- a/tests/test_coreg/test_affine.py +++ b/tests/test_coreg/test_affine.py @@ -2,6 +2,7 @@ from __future__ import annotations +import os.path import warnings import geopandas as gpd @@ -11,7 +12,6 @@ import pytransform3d import rasterio as rio from geoutils import Raster, Vector -from geoutils._typing import NDArrayNum from geoutils.raster import RasterType from geoutils.raster.geotransformations import _translate from scipy.ndimage import binary_dilation @@ -42,53 +42,6 @@ def load_examples(crop: bool = True) -> tuple[RasterType, RasterType, Vector]: return reference_dem, to_be_aligned_dem, glacier_mask -def gdal_reproject_horizontal_shift_samecrs(filepath_example: str, xoff: float, yoff: float) -> NDArrayNum: - """ - Reproject horizontal shift in same CRS with GDAL for testing purposes. - - :param filepath_example: Path to raster file. - :param xoff: X shift in georeferenced unit. - :param yoff: Y shift in georeferenced unit. - - :return: Reprojected shift array in the same CRS. - """ - - from osgeo import gdal, gdalconst - - # Open source raster from file - src = gdal.Open(filepath_example, gdalconst.GA_ReadOnly) - - # Create output raster in memory - driver = "MEM" - method = gdal.GRA_Bilinear - drv = gdal.GetDriverByName(driver) - dest = drv.Create("", src.RasterXSize, src.RasterYSize, 1, gdal.GDT_Float32) - proj = src.GetProjection() - ndv = src.GetRasterBand(1).GetNoDataValue() - dest.SetProjection(proj) - - # Shift the horizontally shifted geotransform - gt = src.GetGeoTransform() - gtl = list(gt) - gtl[0] += xoff - gtl[3] += yoff - dest.SetGeoTransform(tuple(gtl)) - - # Copy the raster metadata of the source to dest - dest.SetMetadata(src.GetMetadata()) - dest.GetRasterBand(1).SetNoDataValue(ndv) - dest.GetRasterBand(1).Fill(ndv) - - # Reproject with resampling - gdal.ReprojectImage(src, dest, proj, proj, method) - - # Extract reprojected array - array = dest.GetRasterBand(1).ReadAsArray().astype("float32") - array[array == ndv] = np.nan - - return array - - class TestAffineCoreg: ref, tba, outlines = load_examples() # Load example reference, to-be-aligned and mask. @@ -121,7 +74,7 @@ class TestAffineCoreg: "xoff_yoff", [(ref.res[0], ref.res[1]), (10 * ref.res[0], 10 * ref.res[1]), (-1.2 * ref.res[0], -1.2 * ref.res[1])], ) # type: ignore - def test_reproject_horizontal_shift_samecrs__gdal(self, xoff_yoff: tuple[float, float]) -> None: + def test_reproject_horizontal_shift_samecrs__gdal(self, xoff_yoff: tuple[float, float], get_test_data_path) -> None: """Check that the same-CRS reprojection based on SciPy (replacing Rasterio due to subpixel errors) is accurate by comparing to GDAL.""" @@ -135,7 +88,8 @@ def test_reproject_horizontal_shift_samecrs__gdal(self, xoff_yoff: tuple[float, ) # Reproject with GDAL - output2 = gdal_reproject_horizontal_shift_samecrs(filepath_example=ref.filename, xoff=xoff, yoff=yoff) + path_output2 = get_test_data_path(os.path.join("gdal", f"shifted_reprojected_xoff{xoff}_yoff{yoff}.tif")) + output2 = Raster(path_output2).data.data # Reproject and NaN propagation is exactly the same for shifts that are a multiple of pixel resolution if xoff % ref.res[0] == 0 and yoff % ref.res[1] == 0: diff --git a/tests/test_terrain.py b/tests/test_terrain.py index 933ec224..52ce4088 100644 --- a/tests/test_terrain.py +++ b/tests/test_terrain.py @@ -1,8 +1,7 @@ from __future__ import annotations -import os +import os.path import re -import tempfile import warnings import geoutils as gu @@ -10,50 +9,12 @@ import pytest import xdem -from xdem._typing import MArrayf xdem.examples.download_longyearbyen_examples() PLOT = True -def run_gdaldem(filepath: str, processing: str, options: str | None = None) -> MArrayf: - """Run GDAL's DEMProcessing and return the read numpy array.""" - # Rasterio strongly recommends against importing gdal along rio, so this is done here instead. - from osgeo import gdal - - gdal.UseExceptions() - - # Converting string into gdal processing options here to avoid import gdal outside this function: - # Riley or Wilson for Terrain Ruggedness, and Zevenberg or Horn for slope, aspect and hillshade - gdal_option_conversion = { - "Riley": gdal.DEMProcessingOptions(alg="Riley"), - "Wilson": gdal.DEMProcessingOptions(alg="Wilson"), - "Zevenberg": gdal.DEMProcessingOptions(alg="ZevenbergenThorne"), - "Horn": gdal.DEMProcessingOptions(alg="Horn"), - "hillshade_Zevenberg": gdal.DEMProcessingOptions(azimuth=315, altitude=45, alg="ZevenbergenThorne"), - "hillshade_Horn": gdal.DEMProcessingOptions(azimuth=315, altitude=45, alg="Horn"), - } - - if options is None: - gdal_option = gdal.DEMProcessingOptions(options=None) - else: - gdal_option = gdal_option_conversion[options] - - temp_dir = tempfile.TemporaryDirectory() - temp_path = os.path.join(temp_dir.name, "output.tif") - gdal.DEMProcessing( - destName=temp_path, - srcDS=filepath, - processing=processing, - options=gdal_option, - ) - - data = gu.Raster(temp_path).data - temp_dir.cleanup() - return data - - class TestTerrainAttribute: filepath = xdem.examples.get_path("longyearbyen_ref_dem") @@ -76,7 +37,7 @@ class TestTerrainAttribute: "roughness", ], ) # type: ignore - def test_attribute_functions_against_gdaldem(self, attribute: str) -> None: + def test_attribute_functions_against_gdaldem(self, attribute: str, get_test_data_path) -> None: """ Test that all attribute functions give the same results as those of GDALDEM within a small tolerance. @@ -100,30 +61,12 @@ def test_attribute_functions_against_gdaldem(self, attribute: str) -> None: "roughness": lambda dem: xdem.terrain.roughness(dem.data), } - # Writing dictionary options here to avoid importing gdal outside the dedicated function - gdal_processing_attr_option = { - "slope_Horn": ("slope", "Horn"), - "aspect_Horn": ("aspect", "Horn"), - "hillshade_Horn": ("hillshade", "hillshade_Horn"), - "slope_Zevenberg": ("slope", "Zevenberg"), - "aspect_Zevenberg": ("aspect", "Zevenberg"), - "hillshade_Zevenberg": ("hillshade", "hillshade_Zevenberg"), - "tri_Riley": ("TRI", "Riley"), - "tri_Wilson": ("TRI", "Wilson"), - "tpi": ("TPI", None), - "roughness": ("Roughness", None), - } - # Copy the DEM to ensure that the inter-test state is unchanged, and because the mask will be modified. dem = self.dem.copy() # Derive the attribute using both GDAL and xdem attr_xdem = functions[attribute](dem).squeeze() - attr_gdal = run_gdaldem( - self.filepath, - processing=gdal_processing_attr_option[attribute][0], - options=gdal_processing_attr_option[attribute][1], - ) + attr_gdal = gu.Raster(get_test_data_path(os.path.join("gdal", f"{attribute}.tif"))).data # For hillshade, we round into an integer to match GDAL's output if attribute in ["hillshade_Horn", "hillshade_Zevenberg"]: @@ -183,7 +126,7 @@ def test_attribute_functions_against_gdaldem(self, attribute: str) -> None: "attribute", ["slope_Horn", "aspect_Horn", "hillshade_Horn", "curvature", "profile_curvature", "planform_curvature"], ) # type: ignore - def test_attribute_functions_against_richdem(self, attribute: str, get_terrain_attribute_richdem) -> None: + def test_attribute_functions_against_richdem(self, attribute: str, get_test_data_path) -> None: """ Test that all attribute functions give the same results as those of RichDEM within a small tolerance. @@ -200,24 +143,13 @@ def test_attribute_functions_against_richdem(self, attribute: str, get_terrain_a "planform_curvature": lambda dem: xdem.terrain.planform_curvature(dem.data, resolution=dem.res), } - # Functions for RichDEM wrapper methods - functions_richdem = { - "slope_Horn": lambda dem: get_terrain_attribute_richdem(dem, attribute="slope", degrees=True), - "aspect_Horn": lambda dem: get_terrain_attribute_richdem(dem, attribute="aspect", degrees=True), - "hillshade_Horn": lambda dem: get_terrain_attribute_richdem(dem, attribute="hillshade"), - "curvature": lambda dem: get_terrain_attribute_richdem(dem, attribute="curvature"), - "profile_curvature": lambda dem: get_terrain_attribute_richdem(dem, attribute="profile_curvature"), - "planform_curvature": lambda dem: get_terrain_attribute_richdem( - dem, attribute="planform_curvature", degrees=True - ), - } - # Copy the DEM to ensure that the inter-test state is unchanged, and because the mask will be modified. dem = self.dem.copy() # Derive the attribute using both RichDEM and xdem attr_xdem = gu.raster.get_array_and_mask(functions_xdem[attribute](dem))[0].squeeze() - attr_richdem = gu.raster.get_array_and_mask(functions_richdem[attribute](dem))[0].squeeze() + attr_richdem_rst = gu.Raster(get_test_data_path(os.path.join("richdem", f"{attribute}.tif")), load_data=True) + attr_richdem = gu.raster.get_array_and_mask(attr_richdem_rst)[0].squeeze() # We compute the difference and keep only valid values diff = attr_xdem - attr_richdem @@ -260,13 +192,13 @@ def test_attribute_functions_against_richdem(self, attribute: str, get_terrain_a raise exception # Introduce some nans - rng = np.random.default_rng(42) - dem.data.mask = np.zeros_like(dem.data, dtype=bool) - dem.data.mask.ravel()[rng.choice(dem.data.size, 50000, replace=False)] = True + # rng = np.random.default_rng(42) + # dem.data.mask = np.zeros_like(dem.data, dtype=bool) + # dem.data.mask.ravel()[rng.choice(dem.data.size, 50000, replace=False)] = True # Validate that this doesn't raise weird warnings after introducing nans and that mask is preserved - output = functions_richdem[attribute](dem) - assert np.all(dem.data.mask == output.data.mask) + # output = functions_richdem[attribute](dem) + # assert np.all(dem.data.mask == output.data.mask) def test_hillshade_errors(self) -> None: """Validate that the hillshade function raises appropriate errors."""