Apollo updates

sdunesme · sdunesme · commit fd72a3c04e76 · 2024-06-10T17:29:17.000+02:00
diff --git a/.gitignore b/.gitignore
@@ -13,6 +13,7 @@ build
 /site
 __pycache__
 .DS_Store
+tutorials/**/*.local.py
 tutorials/outputs/*
 tutorials/inputs/DEM_5m.tif
 !tutorials/outputs/project.rs.xml
diff --git a/fct/cli/Tiles.py b/fct/cli/Tiles.py
@@ -153,11 +153,11 @@ def CreateTileset(datasource: str = 'bdalti',
     
     # Tileset 1
     
-    minx -= resolution/2
-    miny -= resolution/2
+    minx -= (resolution)
+    miny -= (resolution)
     
-    maxx += resolution/2
-    maxy += resolution/2
+    maxx += (resolution)
+    maxy += (resolution)
     
     gx, gy = np.arange(minx, maxx, resolution), np.arange(miny, maxy, resolution)
 
@@ -188,6 +188,7 @@ def CreateTileset(datasource: str = 'bdalti',
     
     minx -= (resolution/2)
     miny -= (resolution/2)
+
     maxx += (resolution/2)
     maxy += (resolution/2)
     
diff --git a/fct/corridor/MedialAxis2.py b/fct/corridor/MedialAxis2.py
@@ -24,21 +24,21 @@
 )
 from shapely.ops import linemerge
 
-from ..config import (
+from fct.config import (
     config,
     LiteralParameter,
     DatasetParameter
 )
-from ..corridor.ValleyBottomFeatures import (
+from fct.corridor.ValleyBottomFeatures import (
     MASK_EXTERIOR,
     MASK_FLOOPLAIN_RELIEF,
     MASK_VALLEY_BOTTOM,
     MASK_TERRACE,
     MASK_SLOPE,
     MASK_HOLE
 )
-from ..simplify import simplify
-from .. import (
+from fct.simplify import simplify
+from fct import (
     transform,
     speedup
 )
@@ -334,8 +334,32 @@ def open_medialaxis_sink(output, crs, driver='ESRI Shapefile'):
                 continue
 
             # lines = list(medial_segments(voronoi, groups[axis]))
-            lines, midpoints, distances = zip(*medial_segments(voronoi, groups[axis]))
             
+            try:
+                tmp_lines, tmp_midpoints, tmp_distances = zip(*medial_segments(voronoi, groups[axis]))
+            except:
+                click.echo(f'Error with axis {axis}. skipping')
+                continue
+            
+            lines = list()
+            midpoints = list()
+            distances = list()
+
+            # 0. remove lines between two axes
+            with rio.open(params.nearest.filename()) as ds:
+                for line, midpoint, distance in zip(tmp_lines, tmp_midpoints, tmp_distances):
+                    line_coords = list(line.coords)
+                    nearest_values = np.array(list(map(int, ds.sample(line_coords, 1))))
+
+                    if np.all(nearest_values == axis) or line.length < (2*params.swath_length):
+                        lines.append(line)
+                        midpoints.append(midpoint)
+                        distances.append(distance)
+
+            lines = tuple(lines)
+            midpoints = tuple(midpoints)
+            distances = tuple(distances)
+                
             # 1. compute medial axis
 
             medialaxis = linemerge(lines)
diff --git a/fct/corridor/ValleyBottomLandcover.py b/fct/corridor/ValleyBottomLandcover.py
@@ -16,6 +16,7 @@
 import os
 from multiprocessing import Pool
 import click
+import cv2
 import rasterio as rio
 from ..config import DatasetParameter
 from ..cli import starcall
@@ -89,7 +90,14 @@ def ValleyBottomLandcoverTile(row, col, params):
     with rio.open(mask_tile) as ds:
 
         mask = ds.read(1)
-        data[mask == ds.nodata] = nodata
+
+        if mask.shape != data.shape:
+            mask = cv2.resize(mask, dsize=(data.shape[1], data.shape[0]), interpolation=cv2.INTER_NEAREST)
+
+        if mask.shape == data.shape:
+            data[mask == ds.nodata] = nodata
+        else:
+            print(f'unable to fit mask and {raster_tile}')
 
     with rio.open(output, 'w', **profile) as dst:
         dst.write(data, 1)
diff --git a/fct/drainage/Accumulate.py b/fct/drainage/Accumulate.py
@@ -327,7 +327,7 @@ def CreateOutletsGraph(params, exterior='exterior-inlets', tileset='default'):
                 # if exterior_flow and os.path.exists(exterior_flow):
                 if not params.exterior_flow.none:
 
-                    exterior_flow = params.exterior_flow.filename(tileset)
+                    exterior_flow = params.exterior_flow.filename()
 
                     with fiona.open(exterior_flow) as fs:
                         for feature in fs:
@@ -337,6 +337,8 @@ def CreateOutletsGraph(params, exterior='exterior-inlets', tileset='default'):
                             if not all([loci >= 0, loci < height, locj >= 0, locj < width]):
                                 continue
 
+                            click.echo(f"exterior-inlet with {feature['properties']['AREAKM2']}km2 in tile {row}_{col}")
+
                             # connect exterior->inlet
 
                             i, j = dem.index(*feature['geometry']['coordinates'])
diff --git a/fct/drainage/Burn.py b/fct/drainage/Burn.py
@@ -50,14 +50,14 @@
 #                 for feature in features:
 #                     fst.write(feature)
 
-def BurnTile(params, row, col, burn_delta=0.0):
+def BurnTile(params, row, col, tileset='default'):
     """
     DOCME
     """
 
-    elevation_raster = params.elevations.tilename(row=row, col=col)
+    elevation_raster = params.elevations.tilename(row=row, col=col, tileset=tileset)
     # config.tileset().tilename(dataset, row=row, col=col)
-    hydrography = params.hydrography.tilename(row=row, col=col)
+    hydrography = params.hydrography.tilename(row=row, col=col, tileset=tileset)
     # config.tileset().tilename('stream-network-draped', row=row, col=col)
 
     with rio.open(elevation_raster) as ds:
@@ -76,7 +76,7 @@ def BurnTile(params, row, col, burn_delta=0.0):
                     for a, b in zip(geom[:-1], geom[1:]):
                         for px, py, z in rasterize_linestringz(a, b):
                             if all([py >= 0, py < height, px >= 0, px < width, not np.isinf(z)]):
-                                elevations[py, px] = z - burn_delta
+                                elevations[py, px] = z - params.offset
         else:
 
             click.secho('File not found: %s' % hydrography, fg='yellow')
diff --git a/fct/drainage/DepressionFill.py b/fct/drainage/DepressionFill.py
@@ -133,7 +133,7 @@ def step(msg):
         if offset < 0:
             elevations = ds.read(1)
         else:
-            elevations = BurnTile(params, row, col, offset)
+            elevations = BurnTile(params, row, col, tileset=tileset)
 
     step('Label flats')
 
diff --git a/fct/drainage/Drape.py b/fct/drainage/Drape.py
@@ -117,7 +117,8 @@ def DrapeNetworkAndAdjustElevations(params):
                         # Clamp z to upstream elevation
                         coordinates[k, 2] = zmin
 
-                    fst.write(feature)
+                    if not np.isinf(feature['geometry']['coordinates']).any():
+                        fst.write(feature)
                     progress.update(1)
 
                     nodez[b] = zmin
@@ -126,9 +127,9 @@ def DrapeNetworkAndAdjustElevations(params):
                     if indegree[node] == 0:
                         queue.append(node)
 
-def SplitStreamNetworkIntoTiles(params):
+def SplitStreamNetworkIntoTiles(params, tileset='default'):
 
-    tileset = config.tileset()
+    tileset_data = config.tileset(tileset)
     networkfile = params.draped.filename(tileset=None)
     # config.filename('stream-network-draped') # filename ok
 
@@ -146,10 +147,10 @@ def SplitStreamNetworkIntoTiles(params):
 
         options = dict(driver=fs.driver, crs=fs.crs, schema=schema)
 
-        with click.progressbar(tileset.tiles(), length=len(tileset)) as iterator:
+        with click.progressbar(tileset_data.tiles(), length=len(tileset_data)) as iterator:
             for tile in iterator:
 
-                output = params.draped.tilename(row=tile.row, col=tile.col)
+                output = params.draped.tilename(row=tile.row, col=tile.col, tileset=tileset)
                 # config.tileset().tilename('stream-network-draped', row=tile.row, col=tile.col)
                 with fiona.open(output, 'w', **options) as dst:
 
@@ -161,7 +162,7 @@ def SplitStreamNetworkIntoTiles(params):
 
                         if intersection.geometryType() == 'LineString':
 
-                            props = feature['properties'].copy()
+                            props = feature['properties']
                             props.update(ROW=tile.row, COL=tile.col)
                             dst.write({
                                 'geometry': intersection.__geo_interface__,
@@ -173,7 +174,7 @@ def SplitStreamNetworkIntoTiles(params):
                             for geom in intersection.geoms:
                                 if geom.geometryType() == 'LineString':
 
-                                    props = feature['properties'].copy()
+                                    props = feature['properties']
                                     props.update(ROW=tile.row, COL=tile.col)
                                     dst.write({
                                         'geometry': geom.__geo_interface__,
diff --git a/fct/drainage/FixNoFlow.py b/fct/drainage/FixNoFlow.py
@@ -416,7 +416,12 @@ def write_tile(row, col): #pylint:disable=unused-argument
             x, y = ds.xy(i, j)
             # row, col = xy2tile(x, y)
             row, col = config.tileset(tileset1).index(x, y)
-            read_tile(row, col)
+
+            try:
+                read_tile(row, col)
+            except rio.RasterioIOError:
+                raise ValueError('No match for (%f, %f)' % (x0, y0))
+
             i, j = ds.index(x, y)
 
     # step3. walk downstream on drainage 2
@@ -440,7 +445,12 @@ def write_tile(row, col): #pylint:disable=unused-argument
             x, y = ds.xy(i, j)
             # row, col = xy2tile(x, y)
             row, col = config.tileset(tileset1).index(x, y)
-            read_tile(row, col)
+
+            try:
+                read_tile(row, col)
+            except rio.RasterioIOError:
+                raise ValueError('No match for (%f, %f)' % (x0, y0))
+            
             i, j = ds.index(x, y)
 
     write_tile(row, col)
diff --git a/fct/drainage/IdentifyNetworkNodes.py b/fct/drainage/IdentifyNetworkNodes.py
@@ -184,6 +184,9 @@ def JoinSourcesAttributes(params, tileset='default'):
         schema = fs.schema
         crs = fiona.crs.from_epsg(config.srid)
         options = dict(driver=driver, crs=crs, schema=schema)
+
+        if not 'GID' in schema['properties'].keys():
+            schema['properties']['GID'] = 'int:10'
         
         with fiona.open(params.sources_identified.filename(tileset=tileset), 'w', **options) as dst:
             
diff --git a/fct/drainage/JoinNetworkAttributes.py b/fct/drainage/JoinNetworkAttributes.py
@@ -96,7 +96,7 @@ def JoinNetworkAttributes(
 
                 properties = feature['properties']
                 node = properties['GID']
-                axis = properties['AXIS']
+                axis = int(properties['AXIS'])
                 axis_increment = max(axis, axis_increment)
                 sources[node] = properties
 
@@ -442,4 +442,7 @@ def get_nodeb(a, axis):
                         
                     feature['geometry']['coordinates'] = coords
                     
-                    dst.write(feature)
+                    try:
+                        dst.write(feature)
+                    except TypeError as e:
+                        click.secho(f"Error for axis {axis}: {e}", fg = 'yellow')
diff --git a/fct/measure/SimplifySwathPolygons2.py b/fct/measure/SimplifySwathPolygons2.py
@@ -126,7 +126,7 @@ def SimplifySwathPolygons(params):
             for geometry, fid in simplified:
 
                 feature = fs.get(fid)
-                feature.update(geometry=geometry.__geo_interface__)
+                feature.update(geometry=fiona.Geometry.from_dict(geometry.__geo_interface__))
                 dst.write(feature)
 
 # def SmoothSwathPolygons(
diff --git a/fct/measure/SwathMeasurement.py b/fct/measure/SwathMeasurement.py
@@ -32,19 +32,19 @@
 import fiona.crs
 from shapely.geometry import shape, Polygon
 
-from ..config import (
+from fct.config import (
     config,
     LiteralParameter,
     DatasetParameter
 )
-# from ..tileio import ReadRasterTile
-from ..tileio import as_window
-from ..rasterize import rasterize_linestring, rasterize_linestringz
-from .. import transform as fct
-from .. import speedup
-from .. import terrain_analysis as ta
-from ..cli import starcall
-from ..metadata import set_metadata
+# from fct.tileio import ReadRasterTile
+from fct.tileio import as_window
+from fct.rasterize import rasterize_linestring, rasterize_linestringz
+from fct import transform as fct
+from fct import speedup
+from fct import terrain_analysis as ta
+from fct.cli import starcall
+from fct.metadata import set_metadata
 
 def nearest_value_and_distance(refpixels, domain, nodata):
     """
@@ -306,7 +306,7 @@ def accept(i, j):
 
         mmin = math.floor(mmin / mdelta) * mdelta
         mmax = math.ceil(mmax / mdelta) * mdelta
-        breaks = np.arange(mmin, mmax + mdelta, mdelta)
+        breaks = np.arange(mmin, mmax + 1.5*mdelta, mdelta)
 
         # click.echo('Calculate Measure & Distance Raster')
 
diff --git a/fct/measure/SwathPolygons.py b/fct/measure/SwathPolygons.py
@@ -9,6 +9,7 @@
 from multiprocessing import Pool
 import logging
 
+import fiona.errors
 import numpy as np
 import click
 
@@ -18,20 +19,20 @@
 import fiona.crs
 from shapely.geometry import shape, Polygon
 
-from ..config import (
+from fct.config import (
     config,
     LiteralParameter,
     DatasetParameter
 )
-# from ..tileio import ReadRasterTile
-from ..tileio import as_window
+# from fct.tileio import ReadRasterTile
+from fct.tileio import as_window
 
-from .. import transform as fct
-from .. import speedup
-from ..cli import starcall
+from fct import transform as fct
+from fct import speedup
+from fct.cli import starcall
 
-from ..corridor.SwathDrainage import create_interpolate_drainage_fun
-from ..corridor.ValleyBottomFeatures import (
+from fct.corridor.SwathDrainage import create_interpolate_drainage_fun
+from fct.corridor.ValleyBottomFeatures import (
     MASK_FLOOPLAIN_RELIEF,
     MASK_VALLEY_BOTTOM,
     MASK_HOLE
@@ -372,7 +373,7 @@ def arguments():
         'geometry': 'Polygon',
         'properties': [
             ('GID', 'int'),
-            ('AXIS', 'int:4'),
+            ('AXIS', 'int64:10'),
             ('VALUE', 'int:4'),
             # ('ROW', 'int:3'),
             # ('COL', 'int:3'),
@@ -434,8 +435,11 @@ def arguments():
                                 'DRAINAGE': float(drainage_area)
                             }
                         }
-
-                        dst.write(feature)
+                    
+                        try:
+                            dst.write(feature)
+                        except:
+                            click.secho('Invalid swath at M %d for axis %d' % (measure, axis), fg='red')
 
                         for ring in geom.interiors:
 
diff --git a/fct/profiles/SwathProfile.py b/fct/profiles/SwathProfile.py
diff --git a/tutorials/03-metrics/03-metrics.py b/tutorials/03-metrics/03-metrics.py
