Reshape dataset when landcover resolution != dem resolution and smooth medialaxis

Author: sdunesme

fct/__init__.py

@@ -11,4 +11,4 @@
 ***************************************************************************
 """
 
-__version__ = '1.1.0'
+__version__ = '1.2.0'

fct/cli/Tiles.py

@@ -18,13 +18,15 @@
 import click
 
 import rasterio as rio
+from rasterio.windows import from_bounds
 import numpy as np
 import fiona
 import fiona.crs
 
-from ..config import config
-from ..tileio import as_window
-from ..cli import starcall
+from fct.config import config
+from fct.config.descriptors import DatasetResolver
+from fct.tileio import as_window
+from fct.cli import starcall
 
 def ExtractTile(datasource, dataset, tile, tileset, overwrite=False):
 
@@ -35,10 +37,13 @@ def ExtractTile(datasource, dataset, tile, tileset, overwrite=False):
         return
 
     with rio.open(raster) as ds:
+        #window = as_window(tile.bounds, ds.transform)
+        window = from_bounds(*tile.bounds, ds.transform)
 
-        window = as_window(tile.bounds, ds.transform)
         data = ds.read(1, window=window, boundless=True, fill_value=ds.nodata)
-        transform = ds.transform * ds.transform.translation(window.col_off, window.row_off)
+        
+        # transform = ds.transform * ds.transform.translation(window.col_off, window.row_off)
+        transform = ds.window_transform(window)
 
         profile = ds.profile.copy()
         profile.update(

fct/corridor/ContinuityAnalysisMax.py

@@ -19,24 +19,24 @@
 import logging
 from multiprocessing import Pool
 import numpy as np
+import cv2
 
 import click
 import rasterio as rio
 from rasterio.windows import Window
 
-from .. import speedup
-from ..cli import starcall
-from ..config import (
+from fct import speedup
+from fct.cli import starcall
+from fct.config import (
     config,
     LiteralParameter,
     DatasetParameter
 )
-from .. import transform as fct
-from ..tileio import (
+from fct import transform as fct
+from fct.tileio import (
     PadRaster,
     border
 )
-from .ValleyBottomFeatures import MASK_EXTERIOR
 
 class Parameters:
     """
@@ -154,10 +154,16 @@ def BoundlessRaster(row, col, dataset):
 
     else:
 
+        if nearest_distance.shape != landcover.shape:
+            nearest_distance = cv2.resize(nearest_distance, landcover.shape, interpolation=cv2.INTER_NEAREST)
+
         out = np.full_like(landcover, nodata)
         distance = np.zeros_like(nearest_distance, dtype='float32')
         state = np.uint8(np.abs(nearest_distance) < 1)
 
+        if state.shape != landcover.shape:
+            click.secho(f'Unable to fit tile {row}_{col}', fg='red')
+        
     state[landcover == nodata] = 255
 
     if seeds:

fct/corridor/MedialAxis2.py

@@ -19,10 +19,10 @@
 
 from shapely.geometry import (
     LineString,
-    MultiLineString,
     shape
 )
 from shapely.ops import linemerge
+from shapelysmooth import chaikin_smooth
 
 from fct.config import (
     config,
@@ -289,6 +289,24 @@ def FixOrientation(medialaxis, params):
 
     return medialaxis
 
+
+def OrderAndMergeLines(geoms, params):
+    """
+    Order list of LineString from upstream to downstream
+    """
+    measure_raster = params.measure.filename()
+
+    coords = [list(l.coords) for l in geoms]
+
+    with rio.open(measure_raster) as ds:
+        heights = [(np.median(list(ds.sample(c))).item(), c) for c in coords]
+        heights = sorted(heights, key=lambda x: x[0], reverse=True)
+        
+        coords = [l[1] for l in heights if l[0] != ds.nodata and l[0] > 0]
+        medialaxis = LineString([i for sublist in coords for i in sublist])
+
+    return medialaxis
+
 def MedialAxis(params):
     """
     Calculate corridor medial axis
@@ -380,7 +398,13 @@ def open_medialaxis_sink(output, crs, driver='ESRI Shapefile'):
                     continue
 
                 lines = sorted(lines, key=lambda line: -line.length)
-                medialaxis = lines[0]
+                # medialaxis = lines[0]
+                filtered = [l for l in lines if l.length > 20*params.swath_length]
+
+                if not filtered:
+                    medialaxis = lines[0]
+                else:
+                    medialaxis = linemerge(filtered)
 
             else:
 
@@ -399,7 +423,8 @@ def open_medialaxis_sink(output, crs, driver='ESRI Shapefile'):
 
                     geoms.append(FixOrientation(geom, params))
 
-                medialaxis = MultiLineString(geoms)
+                medialaxis = OrderAndMergeLines(geoms, params)
+                # medialaxis = MultiLineString(geoms)
 
             else:
 
@@ -551,21 +576,15 @@ def open_output_sink(filename, options):
 
             widths = np.interp(measures, smoothed.measure, smoothed.values)
 
-            # simplified = LineString(
-            #     smooth_chaikin(
-            #         np.array([
-            #             c for c, weight, width
-            #             in zip(coords, weights, widths)
-            #             if filtr(weight, width)
-            #         ])))
-
             simplified = LineString([
                 coord for coord, weight, width
                 in zip(coords, weights, widths)
                 if filtr(weight, width)
             ])
 
+            smoothed = chaikin_smooth(simplified, iters=2)
+
             # print(len(simplified.coords))
-            sink.send((axis, simplified))
+            sink.send((axis, smoothed))
 
     sink.close()

fct/profiles/SwathProfile.py

@@ -152,7 +152,7 @@ def SwathProfileUnit(
         nearest = ds.read(1, window=window, boundless=True, fill_value=ds.nodata)
 
         if nearest.shape != values.shape:
-            nearest = cv2.resize(nearest.astype('float32'), dsize=(values.shape[1],values.shape[0]), interpolation=cv2.INTER_NEAREST).astype('uint32')
+            nearest = cv2.resize(nearest.astype('float32'), dsize=(values.shape[1],values.shape[0]), interpolation=cv2.INTER_NEAREST)
 
         # nearest_nodata = ds.nodata
 
@@ -183,7 +183,7 @@ def SwathProfileUnit(
         swaths = measure_to_swath_identifier(swaths, params.swath_length)
 
         if swaths.shape != values.shape:
-            swaths = cv2.resize(swaths.astype('float32'), dsize=(values.shape[1],values.shape[0]), interpolation=cv2.INTER_NEAREST).astype('uint32')
+            swaths = cv2.resize(swaths.astype('float32'), dsize=(values.shape[1],values.shape[0]), interpolation=cv2.INTER_NEAREST)
 
     try:
 
@@ -197,7 +197,7 @@ def SwathProfileUnit(
         logger.error('Window error on swath (%d, %f)', axis, measure)
         return None
 
-    mask = (nearest == axis) & (swaths == swath) & (values != nodata)
+    mask = (nearest == float(axis)) & (swaths == float(swath)) & (values != nodata)
 
     if np.sum(mask) == 0:
         return None

fct/profiles/ValleyBottomHeight.py

@@ -4,24 +4,17 @@
 
 from multiprocessing import Pool
 import numpy as np
-from scipy.interpolate import interp1d
-from scipy.linalg import lstsq
 from sklearn.linear_model import HuberRegressor
-from sklearn.preprocessing import StandardScaler
 from sklearn.exceptions import ConvergenceWarning
 
 import click
-import fiona
 import rasterio as rio
-from shapely.geometry import LineString, asShape
 import xarray as xr
 
-from ..cli import starcall
-from ..tileio import as_window
-from ..config import DatasetParameter
-from ..corridor.ValleyBottomFeatures import MASK_VALLEY_BOTTOM
-from ..metadata import set_metadata
-from .. import speedup
+from fct.cli import starcall
+from fct.tileio import as_window
+from fct.config import DatasetParameter
+from fct.corridor.ValleyBottomFeatures import MASK_VALLEY_BOTTOM
 
 class Parameters:
     """
@@ -233,7 +226,13 @@ def values():
                 bounds,
                 params)
 
-    return xr.concat(values(), 'swath', 'all')
+    try:
+        vals = xr.concat(values(), 'swath', 'all')
+    except ValueError:
+        click.secho(f'No data for axis {axis}', fg='red')
+        vals = None
+
+    return vals
 
 def ValleyBottomHeight(swath_bounds, params: Parameters, processes: int = 1, **kwargs):
 
@@ -286,4 +285,5 @@ def arguments():
         with click.progressbar(pooled, length=length) as iterator:
             values = list(iterator)
 
+    values = [v for v in values if v is not None]
     return xr.concat(values, 'swath', 'all')

fct/tileio.py

@@ -10,7 +10,7 @@
 
 import numpy as np
 import rasterio as rio
-from rasterio.windows import Window
+from rasterio.windows import Window, from_bounds
 from rasterio.warp import Resampling
 
 from .config import config
@@ -43,21 +43,29 @@ def border(height, width):
             yield i, j
         offset = 0
 
+# def as_window(bounds, transform):
+#     """
+#     Convert real world bounds (minx, miny, maxx, maxy)
+#     to raster window using defined RasterIO geo-transform 
+#     """
+
+#     minx, miny, maxx, maxy = bounds
+
+#     row_offset, col_offset = fct.index(minx, maxy, transform)
+#     row_end, col_end = fct.index(maxx, miny, transform)
+
+#     height = row_end - row_offset
+#     width = col_end - col_offset
+
+#     return Window(col_offset, row_offset, width, height)
+
 def as_window(bounds, transform):
     """
     Convert real world bounds (minx, miny, maxx, maxy)
     to raster window using defined RasterIO geo-transform 
     """
 
-    minx, miny, maxx, maxy = bounds
-
-    row_offset, col_offset = fct.index(minx, maxy, transform)
-    row_end, col_end = fct.index(maxx, miny, transform)
-
-    height = row_end - row_offset
-    width = col_end - col_offset
-
-    return Window(col_offset, row_offset, width, height)
+    return from_bounds(*bounds, transform)
 
 def grow_window(window, padding):
     """

requirements/fct.txt

@@ -2,9 +2,9 @@ cython
 wheel
 numpy==1.23.5
 pyyaml==6.0
-xarray==2023.1.0
+xarray==2024.6.0
 netcdf4==1.6.2
-scipy==1.10.0
+scipy==1.14.0
 rasterio==1.3.5
 fiona==1.9.0
 shapely==1.7.0
@@ -14,3 +14,4 @@ duecredit==0.9.2
 python-dotenv==0.21.1
 matplotlib==3.6.3
 rastachimp==0.3
+shapelysmooth==0.2.0