Merge branch 'main' into remove_deprecated_coreg

.github/workflows/python-package.yml

@@ -0,0 +1,41 @@
+# This workflow will install Python dependencies, run tests and lint with a single version of Python
+# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions
+
+name: build
+
+on:
+  push:
+    branches: [ main ]
+  pull_request:
+    branches: [ main ]
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: [3.7, 3.8, 3.9]
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v2
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        $CONDA/bin/conda install -c conda-forge mamba -y
+        $CONDA/bin/mamba env update --file environment.yml --name base
+        $CONDA/bin/pip install git+https://github.com/GlacioHack/GeoUtils.git
+        $CONDA/bin/pip install -r dev-requirements.txt
+        $CONDA/bin/pip install .
+    - name: Lint with flake8
+      run: |
+        # stop the build if there are Python syntax errors or undefined names
+        $CONDA/bin/flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
+        $CONDA/bin/flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    - name: Test with pytest
+      run: |
+        $CONDA/bin/pytest -rA

.github/workflows/python-publish.yml

@@ -1,7 +1,7 @@
 # This workflow will upload a Python Package using Twine when a release is created
 # For more information see: https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#publishing-to-package-registries
 
-name: Upload Python Package
+name: pypi
 
 on:
   release:
@@ -21,11 +21,11 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install setuptools wheel twine
+        python -m pip install -U setuptools wheel build twine --user
     - name: Build and publish
       env:
         TWINE_USERNAME: ${{ secrets.PYPI_USERNAME }}
         TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }}
       run: |
-        python setup.py sdist bdist_wheel
+        python -m build
         twine upload dist/*

README.md

@@ -4,6 +4,7 @@ Set of tools to manipulate Digital Elevation Models (DEMs)
 More documentation to come!
 
 [![Documentation Status](https://readthedocs.org/projects/xdem/badge/?version=latest)](https://xdem.readthedocs.io/en/latest/?badge=latest)
+[![build](https://github.com/GlacioHack/xdem/actions/workflows/python-package.yml/badge.svg)](https://github.com/GlacioHack/xdem/actions/workflows/python-package.yml)
 
 
 ## Installation ##

dev-requirements.txt

@@ -9,3 +9,4 @@ pyproj
 tqdm
 git+https://github.com/GlacioHack/GeoUtils.git
 scikit-gstat
+flake8

docs/source/coregistration.rst

@@ -93,7 +93,7 @@ The loop is stopped either when the maximum iteration limit is reached, or when
         dem_2009 = xdem.DEM(xdem.examples.FILEPATHS["longyearbyen_ref_dem"])
         dem_1990 = xdem.DEM(xdem.examples.FILEPATHS["longyearbyen_tba_dem"])
         outlines_1990 = gu.Vector(xdem.examples.FILEPATHS["longyearbyen_glacier_outlines"])
-        inlier_mask = outlines_1990.create_mask(dem_2009) != 255
+        inlier_mask = ~outlines_1990.create_mask(dem_2009)
 
         nuth_kaab = xdem.coreg.NuthKaab()
         nuth_kaab.fit(dem_2009.data, dem_1990.data, transform=dem_2009.transform, inlier_mask=inlier_mask)
@@ -102,8 +102,8 @@ The loop is stopped either when the maximum iteration limit is reached, or when
         ddem_pre = (dem_2009.data - dem_1990.data).filled(np.nan).squeeze()
         ddem_post = (dem_2009.data - dem_coreg).filled(np.nan).squeeze()
 
-        nmad_pre = xdem.spatial_tools.nmad(ddem_pre[inlier_mask])
-        nmad_post = xdem.spatial_tools.nmad(ddem_post[inlier_mask])
+        nmad_pre = xdem.spatial_tools.nmad(ddem_pre[inlier_mask.squeeze()])
+        nmad_post = xdem.spatial_tools.nmad(ddem_post[inlier_mask.squeeze()])
 
         vlim = 20
         plt.figure(figsize=(8, 5))

environment.yml

@@ -2,8 +2,8 @@ name: xdem
 channels:
   - conda-forge
 dependencies:
-  - python=3
-  - proj=7.2
+  - python>=3
+  - proj>=7.2
   - geopandas
   - ipython
   - numba
@@ -20,5 +20,6 @@ dependencies:
   - pip:
     - richdem
     - scikit-gstat
+    - pytransform3d
 #    - "--editable=git+https://github.com/GlacioHack/GeoUtils.git#egg=GeoUtils"
 #    - "--editable=git+https://github.com/GlacioHack/xdem.git#egg=xdem"

tests/test_coreg.py

@@ -16,6 +16,7 @@
 import cv2
 import geoutils as gu
 import numpy as np
+import pytest
 
 with warnings.catch_warnings():
     warnings.simplefilter("ignore")

tests/test_dem.py

@@ -1,21 +1,24 @@
 """
 Test functions for DEM class
 """
+import inspect
 import os
-import pytest
 import warnings
-import pyproj
-import inspect
+
 import geoutils.georaster as gr
 import geoutils.satimg as si
-import xdem
 import numpy as np
+import pyproj
+import pytest
+
+import xdem
 from xdem.dem import DEM
 
 xdem.examples.download_longyearbyen_examples(overwrite=False)
 
 DO_PLOT = False
 
+
 class TestDEM:
 
     def test_init(self):
@@ -27,38 +30,37 @@ def test_init(self):
 
         # from filename
         dem = DEM(fn_img)
-        assert isinstance(dem,DEM)
+        assert isinstance(dem, DEM)
 
         # from DEM
         dem2 = DEM(dem)
-        assert isinstance(dem2,DEM)
+        assert isinstance(dem2, DEM)
 
         # from Raster
         r = gr.Raster(fn_img)
         dem3 = DEM(r)
-        assert isinstance(dem3,DEM)
+        assert isinstance(dem3, DEM)
 
         # from SatelliteImage
         img = si.SatelliteImage(fn_img)
         dem4 = DEM(img)
-        assert isinstance(dem4,DEM)
+        assert isinstance(dem4, DEM)
 
-        list_dem = [dem,dem2,dem3,dem4]
+        list_dem = [dem, dem2, dem3, dem4]
 
         attrs = [at for at in gr.default_attrs if at not in ['name', 'dataset_mask', 'driver']]
         all_attrs = attrs + si.satimg_attrs + xdem.dem.dem_attrs
         for attr in all_attrs:
             attrs_per_dem = [idem.__getattribute__(attr) for idem in list_dem]
             assert all(at == attrs_per_dem[0] for at in attrs_per_dem)
 
-        assert np.logical_and.reduce((np.array_equal(dem.data,dem2.data,equal_nan=True),
-                               np.array_equal(dem2.data,dem3.data,equal_nan=True),
-                               np.array_equal(dem3.data,dem4.data,equal_nan=True)))
-
-        assert np.logical_and.reduce((np.all(dem.data.mask==dem2.data.mask),
-                               np.all(dem2.data.mask==dem3.data.mask),
-                               np.all(dem3.data.mask==dem4.data.mask)))
+        assert np.logical_and.reduce((np.array_equal(dem.data, dem2.data, equal_nan=True),
+                                      np.array_equal(dem2.data, dem3.data, equal_nan=True),
+                                      np.array_equal(dem3.data, dem4.data, equal_nan=True)))
 
+        assert np.logical_and.reduce((np.all(dem.data.mask == dem2.data.mask),
+                                      np.all(dem2.data.mask == dem3.data.mask),
+                                      np.all(dem3.data.mask == dem4.data.mask)))
 
     def test_copy(self):
         """
@@ -114,15 +116,15 @@ def test_set_vref(self):
         assert img.vref == 'EGM96'
         assert img.vref_grid == 'us_nga_egm96_15.tif'
         # grid should have priority over name and parse the right vref name
-        img.set_vref(vref_name='WGS84',vref_grid='us_nga_egm96_15.tif')
+        img.set_vref(vref_name='WGS84', vref_grid='us_nga_egm96_15.tif')
         assert img.vref == 'EGM96'
 
         # check for EGM08
         img.set_vref(vref_name='EGM08')
         assert img.vref == 'EGM08'
         assert img.vref_grid == 'us_nga_egm08_25.tif'
         # grid should have priority over name and parse the right vref name
-        img.set_vref(vref_name='best ref in the entire world, or any string',vref_grid='us_nga_egm08_25.tif')
+        img.set_vref(vref_name='best ref in the entire world, or any string', vref_grid='us_nga_egm08_25.tif')
         assert img.vref == 'EGM08'
 
         # check that other existing grids are well detected in the pyproj.datadir
@@ -132,75 +134,74 @@ def test_set_vref(self):
         with pytest.raises(ValueError):
             img.set_vref(vref_grid='the best grid in the entire world, or any non-existing string')
 
-
     def test_to_vref(self):
 
-        #first, some points to test the transform
+        # first, some points to test the transform
 
         # Chile
         lat = 43.70012234
         lng = -79.41629234
         z = 100
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", module="pyproj")
-            #init is deprecated by
-            ellipsoid=pyproj.Proj(init="EPSG:4326") #WGS84 datum ellipsoid height
-            geoid=pyproj.Proj(init="EPSG:4326", geoidgrids='us_nga_egm96_15.tif') #EGM96 geoid in Chile, we expect ~30 m difference
+            # init is deprecated by
+            ellipsoid = pyproj.Proj(init="EPSG:4326")  # WGS84 datum ellipsoid height
+            # EGM96 geoid in Chile, we expect ~30 m difference
+            geoid = pyproj.Proj(init="EPSG:4326", geoidgrids='us_nga_egm96_15.tif')
         transformer = pyproj.Transformer.from_proj(ellipsoid, geoid)
-        z_out = transformer.transform(lng,lat,z)[2]
+        z_out = transformer.transform(lng, lat, z)[2]
 
-        #check final elevation is finite, higher than ellipsoid with less than 40 m difference (typical geoid in Chile)
-        assert np.logical_and.reduce((np.isfinite(z_out),np.greater(z_out,z),np.less(np.abs(z_out-z),40)))
+        # check final elevation is finite, higher than ellipsoid with less than 40 m difference (typical geoid in Chile)
+        assert np.logical_and.reduce((np.isfinite(z_out), np.greater(z_out, z), np.less(np.abs(z_out-z), 40)))
 
-
-        #egm2008
+        # egm2008
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", module="pyproj")
-            #init is deprecated by
-            ellipsoid=pyproj.Proj(init="EPSG:4326") #WGS84 datum ellipsoid height
-            geoid=pyproj.Proj(init="EPSG:4326", geoidgrids='us_nga_egm08_25.tif')
+            # init is deprecated by
+            ellipsoid = pyproj.Proj(init="EPSG:4326")  # WGS84 datum ellipsoid height
+            geoid = pyproj.Proj(init="EPSG:4326", geoidgrids='us_nga_egm08_25.tif')
         transformer = pyproj.Transformer.from_proj(ellipsoid, geoid)
-        z_out = transformer.transform(lng,lat,z)[2]
+        z_out = transformer.transform(lng, lat, z)[2]
 
-        #check final elevation is finite, higher than ellipsoid with less than 40 m difference (typical geoid in Chile)
-        assert np.logical_and.reduce((np.isfinite(z_out),np.greater(z_out,z),np.less(np.abs(z_out-z),40)))
+        # check final elevation is finite, higher than ellipsoid with less than 40 m difference (typical geoid in Chile)
+        assert np.logical_and.reduce((np.isfinite(z_out), np.greater(z_out, z), np.less(np.abs(z_out-z), 40)))
 
-        #geoid2006 for Alaska
+        # geoid2006 for Alaska
         lat = 65
         lng = -140
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", module="pyproj")
-            #init is deprecated by
-            ellipsoid=pyproj.Proj(init="EPSG:4326") #WGS84 datum ellipsoid height
-            geoid=pyproj.Proj(init="EPSG:4326", geoidgrids='us_noaa_geoid06_ak.tif')
+            # init is deprecated by
+            ellipsoid = pyproj.Proj(init="EPSG:4326")  # WGS84 datum ellipsoid height
+            geoid = pyproj.Proj(init="EPSG:4326", geoidgrids='us_noaa_geoid06_ak.tif')
         transformer = pyproj.Transformer.from_proj(ellipsoid, geoid)
-        z_out = transformer.transform(lng,lat,z)[2]
-
-        #check final elevation is finite, lower than ellipsoid with less than 20 m difference (typical geoid in Alaska)
-        assert np.logical_and.reduce((np.isfinite(z_out),np.less(z_out,z),np.less(np.abs(z_out-z),20)))
+        z_out = transformer.transform(lng, lat, z)[2]
 
+        # check final elevation is finite, lower than ellipsoid with less than 20 m difference (typical geoid in Alaska)
+        assert np.logical_and.reduce((np.isfinite(z_out), np.less(z_out, z), np.less(np.abs(z_out-z), 20)))
 
-        #isn1993 for Iceland
+        # isn1993 for Iceland
         lat = 65
         lng = -18
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", module="pyproj")
-            #init is deprecated by
-            ellipsoid=pyproj.Proj(init="EPSG:4326") #WGS84 datum ellipsoid height
-            geoid=pyproj.Proj(init="EPSG:4326", geoidgrids='is_lmi_Icegeoid_ISN93.tif') #Iceland, we expect a ~70m difference
+            # init is deprecated by
+            ellipsoid = pyproj.Proj(init="EPSG:4326")  # WGS84 datum ellipsoid height
+            # Iceland, we expect a ~70m difference
+            geoid = pyproj.Proj(init="EPSG:4326", geoidgrids='is_lmi_Icegeoid_ISN93.tif')
         transformer = pyproj.Transformer.from_proj(ellipsoid, geoid)
-        z_out = transformer.transform(lng,lat,z)[2]
+        z_out = transformer.transform(lng, lat, z)[2]
 
-        #check final elevation is finite, lower than ellipsoid with less than 100 m difference (typical geoid in Iceland)
-        assert np.logical_and.reduce((np.isfinite(z_out),np.less(z_out,z),np.less(np.abs(z_out-z),100)))
+        # check final elevation is finite, lower than ellipsoid with less than 100 m difference (typical geoid in Iceland)
+        assert np.logical_and.reduce((np.isfinite(z_out), np.less(z_out, z), np.less(np.abs(z_out-z), 100)))
 
-        #checking that the function does not run without a reference set
+        # checking that the function does not run without a reference set
         fn_img = xdem.examples.FILEPATHS["longyearbyen_ref_dem"]
         img = DEM(fn_img)
         with pytest.raises(ValueError):
             img.to_vref(vref_name='EGM96')
 
-        #checking that the function properly runs with a reference set
+        # checking that the function properly runs with a reference set
         img.set_vref(vref_name='WGS84')
         mean_ellips = np.nanmean(img.data)
         img.to_vref(vref_name='EGM96')
@@ -209,11 +210,7 @@ def test_to_vref(self):
         assert img.vref == 'EGM96'
         assert img.vref_grid == 'us_nga_egm96_15.tif'
 
-        #check that the geoid is lower than ellipsoid, less than 35 m difference (Svalbard)
-
-        assert np.greater(mean_ellips,mean_geoid_96)
-        assert np.less(np.abs(mean_ellips-mean_geoid_96),35.)
-
-
-
+        # check that the geoid is lower than ellipsoid, less than 35 m difference (Svalbard)
 
+        assert np.greater(mean_ellips, mean_geoid_96)
+        assert np.less(np.abs(mean_ellips-mean_geoid_96), 35.)

tests/test_spatial_tools.py

@@ -52,4 +52,4 @@ def test_merge_rasters():
 
     diff = dem.data - merged_dem.data
 
-    assert np.abs(np.nanmean(diff)) < 0.0001
+    assert np.abs(np.nanmean(diff)) < 0.05