-
Notifications
You must be signed in to change notification settings - Fork 42
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.
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
Remove GDAL and RichDEM dependancy from tests #675
base: main
Are you sure you want to change the base?
Changes from all commits
6c55595
ebe17a6
3baa2fb
1f9e9e7
43b609a
f0a24a7
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -9,7 +9,7 @@ dependencies: | |
- matplotlib=3.* | ||
- pyproj>=3.4,<4 | ||
- rasterio>=1.3,<2 | ||
- scipy=1.* | ||
- scipy<1.15.0 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. same review |
||
- tqdm | ||
- scikit-image=0.* | ||
- scikit-gstat>=1.0.18,<1.1 | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -7,7 +7,7 @@ numpy==1.* | |
matplotlib==3.* | ||
pyproj>=3.4,<4 | ||
rasterio>=1.3,<2 | ||
scipy==1.* | ||
scipy<1.15.0 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. same |
||
tqdm | ||
scikit-image==0.* | ||
scikit-gstat>=1.0.18,<1.1 | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -61,7 +61,6 @@ test = | |
flake8 | ||
pylint | ||
scikit-learn | ||
richdem | ||
doc = | ||
sphinx | ||
sphinx-book-theme | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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" | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I quickly checked if downloading was the best solution, and it seems fine and easy to maintain. Indeed, it might be better to download a tarball via a release. Is that what you were thinking? Also, we obviously need to update this PR once the xdem-data PR is merged and a release is created, so we can have the correct tarball/zip with a fixed version API. It doesn't seem like you're handling a version here. Maybe it's better to use a release, something like: https://github.com/glaciohack/xdem-data/releases/download/v1.0.0/xdem-data.zip There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @duboise-cnes I used the same download system as for importing the example data (Longyearbyen DEMs) here to maintain the same consistency and ease of maintenance. The example data are downloading from the main There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. ok I understand therefore the rationale, you know better the code of xdem, i didn't know that there was the example download code. I would say that using requests is better than urllib for maintenance but these would change the entire also for the examples. And the download function could not be factorized with the example data, some code replication here. Not needed in my opinion. |
||
|
||
@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: | ||
""" | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think that you can use requests et zipfile (or tarball as used equivalent) if a release ? For example from a well known AI ;) not used directly just for inspiration
does it help ? |
||
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: | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. do not needed with requests ? |
||
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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. check the format with a versioned release in xdem-data. if tar ok, otherwise zip equivalent. |
||
|
||
# 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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
why <1.15.0 ? not great to keep that. If really needed, please add an issue to correct it in another PR. This can be difficult with other packages in time. I prefer scipy >=, !=1.15.0 . it allows that next scipy release automatic upgrade and test if bugs are resolved.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
#677