Feature/reimplement nearneighbour

.github/workflows/ci.yml

@@ -20,26 +20,28 @@ jobs:
     steps:
       - uses: actions/checkout@v2
 
-      - name: Setup conda
-        uses: s-weigand/setup-conda@v1
+      - uses: conda-incubator/setup-miniconda@v2
         with:
+          activate-environment: pycascadia
+          environment-file: environment.yml
           python-version: ${{ matrix.python }}
-          conda-channels: conda-forge
-
-      - name: Install conda dependencies
-        run: |
-          conda install pip
-          conda install -c conda-forge pygmt gdal rasterio
+          auto-activate-base: false
 
       - name: Install package
-        run: pip install .[test] flake8
+        shell: bash -l {0}
+        run: pip install .[test]
 
       - name: Lint with flake8
+        shell: bash -l {0}
         run: |
+          conda install flake8
           # stop the build if there are Python syntax errors or undefined names
           flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
           # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
           flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
 
       - name: Test
-        run: pytest
+        shell: bash -l {0}
+        run: |
+          conda install pytest
+          pytest

README.md

@@ -9,27 +9,23 @@ by running the commands below in sequence:
 ```
 git clone https://github.com/UCL/pyCascadia.git
 ```
-2. In your local copy of the GitHub repository, create a conda environment (named "cascadia" here, but you might want to choose a different name) with Python 3.9. Other Python versions might also work, but the code is only actively tested with Python 3.9. This step also installs `pip`, which will be used to install the more straight-forward dependencies.
+2. The proceeding steps assume your current working directory is your local copy of the repository. A mixture of conda and pip is used to install the dependencies (due to the complex dependencies gdal and rasterio). The conda dependencies are listed in [environment.yml](https://github.com/UCL/pyCascadia/blob/main/environment.yml) which will be automatically used by conda when creating a new environment:
 ```
-conda create -n cascadia python=3.9 pip
+conda env create
 ```
-3. Additionally, install `pygmt`, `gdal` and `rasterio` into the "cascadia" environment via `conda-forge`. This step is needed because these packages are not available from the Python Package index (PyPi) (which is where `pip` looks by default) for all operating systems.
+3. This will create a new conda environment named `pycascadia` which should be activated with:
 ```
-conda install -n cascadia pygmt gdal rasterio -c conda-forge
+conda activate pycascadia
 ```
-4. Activate the "cascadia" environment
-```
-conda activate cascadia
-```
-5. Install `pyCascadia`. This will also automatically install the remaining dependencies (a list of which can be found in [setup.cfg](https://github.com/UCL/pyCascadia/blob/main/setup.cfg)).
+4. Install `pyCascadia`. This will also automatically install the remaining dependencies (a list of which can be found in [setup.cfg](https://github.com/UCL/pyCascadia/blob/main/setup.cfg)).
 ```
 pip install .[test]
 ```
-6. Check the installation worked by running the unit tests.
+5. Check the installation worked by running the unit tests:
 ```
 pytest
 ```
-7. Once you've done your work with `pyCascadia`, you may want to deactivate the environment and return to your base python environment. You can do so by running:
+6. Once you've done your work with `pyCascadia`, you may want to deactivate the environment and return to your base python environment. You can do so by running:
 ```
 conda deactivate
 ```

environment.yml

@@ -0,0 +1,11 @@
+name: pycascadia
+channels:
+    - conda-forge
+    - defaults
+dependencies:
+    # Required dependencies
+    - pip
+    - pygmt>=0.3.1
+    - gdal
+    - rasterio
+    - pytest

pycascadia/remove_restore.py

@@ -7,6 +7,16 @@
 """
 
 from pygmt import blockmedian, surface, grdtrack, grdcut, grdfilter
+from pygmt.clib import Session
+from pygmt.helpers import (
+    GMTTempFile,
+    build_arg_string,
+    fmt_docstring,
+    kwargs_to_strings,
+    use_alias,
+)
+import xarray as xr
+
 import os
 import pandas as pd
 import matplotlib.pyplot as plt
@@ -16,6 +26,36 @@
 from pycascadia.grid import Grid
 from pycascadia.utility import region_to_str, min_regions, is_region_valid, read_fnames
 
+@use_alias(
+    I="spacing",
+    R="region",
+    V="verbose",
+)
+@kwargs_to_strings(R="sequence")
+def nearneighbour(data_xyz, **kwargs):
+    """Uses pyGMT's clib to call GMT's nearneighbour command
+
+    Adapted from pyGMT's [blockmedian implementation](https://github.com/GenericMappingTools/pygmt/blob/c0ff7f1add9884305688c2fa15c5f13516b8b960/pygmt/src/blockm.py)
+    """
+    with GMTTempFile(suffix=".csv") as tmpfile:
+        with Session() as lib:
+            file_context = lib.virtualfile_from_data(check_kind="vector", data=data_xyz)
+            with file_context as infile:
+                if "G" not in kwargs.keys():  # if outgrid is unset, output to tempfile
+                    kwargs.update({"G": tmpfile.name})
+                outgrid = kwargs["G"]
+                arg_str = " ".join([infile, build_arg_string(kwargs)])
+                lib.call_module("nearneighbor", arg_str)
+
+        if outgrid == tmpfile.name:  # if user did not set outgrid, return DataArray
+            with xr.open_dataarray(outgrid) as dataarray:
+                result = dataarray.load()
+                _ = result.gmt  # load GMTDataArray accessor information
+        else:
+            result = None  # if user sets an outgrid, return None
+
+    return result
+
 def calc_diff_grid(base_grid, update_grid, diff_threshold=0.0, window_width=None):
     """Calculates difference grid for use in remove-restore"""
     print("Blockmedian update grid")
@@ -38,18 +78,9 @@ def calc_diff_grid(base_grid, update_grid, diff_threshold=0.0, window_width=None
 
     diff[diff.z.abs() < diff_threshold]['z'] = 0.0 # Filter out small differences
 
-    diff_xyz_fname = "diff.xyz"
-    diff_grid_fname = "diff.nc"
-    diff.to_csv(diff_xyz_fname, sep=' ', header=False, index=False)
-
     NODATA_VAL = 9999
 
-    os.system(f'gmt nearneighbor {diff_xyz_fname} -G{diff_grid_fname} -R{region_to_str(base_grid.region)} -I{base_grid.spacing} -S{2*max_spacing} -N4 -E{NODATA_VAL} -V')
-    diff_grid, _, _ = load_source(diff_grid_fname)
-
-    # Cleanup files
-    os.remove(diff_grid_fname)
-    os.remove(diff_xyz_fname)
+    diff_grid = nearneighbour(diff, region=base_grid.region, spacing=base_grid.spacing, S=2*max_spacing, N=4, E=NODATA_VAL, verbose=True)
 
     # Interpolate between nodata and data regions in update grid
     if window_width:

tests/test_remove_restore.py

@@ -0,0 +1,34 @@
+import pytest
+import os
+import pandas as pd
+from xarray.testing import assert_equal
+
+from pycascadia.remove_restore import nearneighbour
+from pycascadia.utility import region_to_str
+from pycascadia.loaders import load_source
+
+def test_nearneighbour():
+    region = [0, 1, 0, 1]
+    spacing = 0.3
+    NODATA_VAL = 9999
+    max_spacing = 0.2
+
+    data_xyz = pd.DataFrame.from_dict({
+        'x': [0.1, 0.4, 0.7],
+        'y': [0.2, 0.3, 0.4],
+        'z': [10, 14, 4],
+    })
+
+    # Create grid via nearneighbour function
+    clib_grid = nearneighbour(data_xyz, region=region, spacing=spacing, S=2*max_spacing, N=4, E=NODATA_VAL, verbose=True)
+
+    # Create grid via manual os call to gmt
+    data_xyz_fname = "data.xyz"
+    data_grid_fname = "data.nc"
+    data_xyz.to_csv(data_xyz_fname, sep=' ', header=False, index=False)
+    os.system(f'gmt nearneighbor {data_xyz_fname} -G{data_grid_fname} -R{region_to_str(region)} -I{spacing} -S{2*max_spacing} -N4 -E{NODATA_VAL} -V')
+    manual_grid, _, _ = load_source(data_grid_fname)
+    os.remove(data_grid_fname)
+    os.remove(data_xyz_fname)
+
+    assert_equal(manual_grid, clib_grid)