Make output directory a keyword argument

topo_descriptors/config/topo_descriptors.conf

@@ -3,7 +3,3 @@ min_elevation: -100
 
 # The number of standard deviations per unit scale
 scale_std: 4
-
-# Paths
-path:
-    output: .

topo_descriptors/helpers.py

@@ -37,7 +37,7 @@ def get_dem_netcdf(path_dem):
     return dem_da.where(dem_da > CFG.min_elevation)
 
 
-def to_netcdf(array, coords, name, crop=None):
+def to_netcdf(array, coords, name, crop=None, outdir="."):
     """Save an array of topographic descriptors in NetCDF. It is first converted
     into a xarray DataArray with the same coordinates as the input DEM DataArray
     and a specified name.
@@ -53,12 +53,15 @@ def to_netcdf(array, coords, name, crop=None):
         The array is cropped to the given extend before being saved. Keys should
         be coordinates labels as in coords and values should be slices of [min,max]
         extend. Default is None.
+    outdir (optional) : string
+        The path to the output directory. Save to working directory by default.
     """
 
     name = str.upper(name)
+    outdir = Path(outdir)
     da = xr.DataArray(array, coords=coords, name=name).sel(crop)
     filename = f"topo_{name}.nc"
-    da.to_dataset().to_netcdf(CFG.path.output / filename)
+    da.to_dataset().to_netcdf(outdir / filename)
 
 
 def scale_to_pixel(scales, dem_da):
@@ -83,7 +86,7 @@ def scale_to_pixel(scales, dem_da):
         Resolution in meters of each DEM grid points in the x and y directions.
     """
     check_dem(dem_da)
-    x_coords, y_coords = dem_da['x'].values, dem_da['y'].values
+    x_coords, y_coords = dem_da["x"].values, dem_da["y"].values
     epsg_code = f"epsg:{int(dem_da.attrs['crs'].split('epsg:')[1])}"
     if epsg_code == "epsg:4326":
         logger.warning(
@@ -176,9 +179,9 @@ def check_dem(dem):
         raise ValueError("dem must be a xr.DataArray")
     if dem.ndim != 2:
         raise ValueError("dem must be a two-dimensional array")
-    if dem.dims != ('y', 'x'):
+    if dem.dims != ("y", "x"):
         raise ValueError("dem dimensions must be ('y', 'x')")
-    if not 'crs' in dem.attrs:
+    if not "crs" in dem.attrs:
         raise KeyError("missing 'crs' attribute in dem")
-    if not 'epsg:' in dem.attrs['crs'].lower():
-        raise ValueError("missing 'epsg:' key in the 'crs' attribute")
+    if not "epsg:" in dem.attrs["crs"].lower():
+        raise ValueError("missing 'epsg:' key in the 'crs' attribute")

topo_descriptors/topo.py

@@ -1,4 +1,5 @@
 import logging
+from pathlib import Path
 
 import numpy as np
 import numpy.ma as ma
@@ -11,7 +12,7 @@
 logger = logging.getLogger(__name__)
 
 
-def compute_tpi(dem_da, scales, smth_factors=None, ind_nans=[], crop=None):
+def compute_tpi(dem_da, scales, smth_factors=None, ind_nans=[], crop=None, outdir="."):
     """Wrapper to 'tpi' function to launch computations for all scales and save
     outputs as netCDF files.
 
@@ -35,6 +36,8 @@ def compute_tpi(dem_da, scales, smth_factors=None, ind_nans=[], crop=None):
         If specified the outputs are cropped to the given extend. Keys should be
         the coordinates labels of dem_da and values should be slices of [min,max]
         extend. Default is None.
+    outdir (optional) : string
+        The path to the output directory. Save to working directory by default.
 
     See also
     --------
@@ -60,7 +63,7 @@ def compute_tpi(dem_da, scales, smth_factors=None, ind_nans=[], crop=None):
         array = tpi(dem=dem_da.values, size=scale_pxl, sigma=sigmas[idx])
 
         array[ind_nans] = np.nan
-        hlp.to_netcdf(array, dem_da.coords, name, crop)
+        hlp.to_netcdf(array, dem_da.coords, name, crop, outdir)
         del array
 
 
@@ -138,6 +141,7 @@ def compute_valley_ridge(
     smth_factors=None,
     ind_nans=[],
     crop=None,
+    outdir=".",
 ):
     """Wrapper to 'valley_ridge' function to launch computations for all scales
     and save outputs as netCDF files.
@@ -168,6 +172,8 @@ def compute_valley_ridge(
         If specified the outputs are cropped to the given extend. Keys should be
         the coordinates labels of dem_da and values should be slices of [min,max]
         extend. Default is None.
+    outdir (optional) : string
+        The path to the output directory. Save to working directory by default.
 
     See also
     --------
@@ -202,21 +208,24 @@ def compute_valley_ridge(
                 names,
                 ind_nans,
                 crop,
+                outdir,
             ),
         )
 
     pool.close()
     pool.join()
 
 
-def _valley_ridge_wrap(dem_da, size, mode, flat_list, sigma, names, ind_nans, crop):
+def _valley_ridge_wrap(
+    dem_da, size, mode, flat_list, sigma, names, ind_nans, crop, outdir
+):
     """Wrapper to valley_ridge and hlp.to_netcdf functions to ease the parallelization
     of the different scales"""
 
     arrays = valley_ridge(dem_da.values, size, mode, flat_list, sigma)
     for array, name in zip(arrays, names):
         array[ind_nans] = np.nan
-        hlp.to_netcdf(array, dem_da.coords, name, crop)
+        hlp.to_netcdf(array, dem_da.coords, name, crop, outdir)
 
 
 @hlp.timer
@@ -365,7 +374,7 @@ def _rotate_kernels(kernel, angle):
     return ma.MaskedArray.filled(kernels_rot, 0).astype(np.float32)
 
 
-def compute_gradient(dem_da, scales, sig_ratios=1, ind_nans=[], crop=None):
+def compute_gradient(dem_da, scales, sig_ratios=1, ind_nans=[], crop=None, outdir="."):
     """Wrapper to 'gradient' function to launch computations for all scales
     and save outputs as netCDF files.
 
@@ -388,6 +397,8 @@ def compute_gradient(dem_da, scales, sig_ratios=1, ind_nans=[], crop=None):
         If specified the outputs are cropped to the given extend. Keys should be
         the coordinates labels of dem_da and values should be slices of [min,max]
         extend. Default is None.
+    outdir (optional) : string
+        The path to the output directory. Save to working directory by default.
 
     See also
     --------
@@ -419,7 +430,7 @@ def compute_gradient(dem_da, scales, sig_ratios=1, ind_nans=[], crop=None):
 
         for array, name in zip(arrays, names):
             array[ind_nans] = np.nan
-            hlp.to_netcdf(array, dem_da.coords, name, crop)
+            hlp.to_netcdf(array, dem_da.coords, name, crop, outdir)
 
         del arrays