AHN from files

sample/geslotenheid/ahn_opgaande_elemenden/ahn_opgaande.py

@@ -8,6 +8,7 @@ try:
     from ahn_read import get_opgaande_elementen, verify_matching_rasters
     from hoogtes_dbf import build_heights_table, heights_table_to_file
     from flt_header import replace_hdr
+    from output_profile import build_output_profile, get_ahn_files
 except ModuleNotFoundError:
     from sample.geslotenheid.ahn_opgaande_elemenden.ahn_read import get_opgaande_elementen, verify_matching_rasters
     from sample.geslotenheid.ahn_opgaande_elemenden.read_write_windows import get_rw_windows
@@ -27,44 +28,66 @@ def build_opgaande_elementenkaart(**kwargs):
 
     """
 
-    dsm_source = kwargs['dsm_src']
-    dtm_source = kwargs['dtm_src']
-
-    dsm = rio.open(dsm_source)
-    dtm = rio.open(dtm_source)
-
-    verify_matching_rasters(dsm, dtm)
-
-    windows = get_rw_windows(target_profile=dsm.profile,
-                             sample=kwargs['sample'],
-                             window_width=kwargs['w_width'],
-                             window_height=kwargs['w_height'])
-
-    output_profile = dsm.profile
-    output_profile.update(dtype=np.uint16, nodata=0)
-    output_profile.update(driver={'tif': "GTiff", 'flt': 'EHdr'}[kwargs['of']])
-    extension = kwargs['of']
+    output_profile, data_source = build_output_profile(**kwargs)
     u_values = set()
 
-    with rio.open(os.path.join(kwargs['out_dir'], f'{kwargs["out_name"]}.{extension}'), 'w', **output_profile) as dest:
+    with rio.open(os.path.join(kwargs['out_dir'], f'{kwargs["out_name"]}.{kwargs["of"]}'), 'w', **output_profile) as dest:
 
         dest.update_tags(**kwargs)
 
-        for i, w in enumerate(windows, start=1):
-            print(f'processing window {i}/{len(windows)}', end="\r")
-            opg = get_opgaande_elementen(
-                dsm=dsm,
-                dtm=dtm,
-                read_window=w.rio_window,
-                minimum_difference=kwargs['min_diff'],
-                uniform_height=kwargs['u_height'],
-                maaiveld_value=kwargs['mv_value'],
-            )
-            dest.write(opg, 1,  window=w.rio_window)
-
-            u_values.update(np.unique(opg[opg != kwargs['mv_value']]))
-
-    heights_table_to_file(os.path.join(kwargs['out_dir'],
+        if data_source == 'file':
+            windows = get_rw_windows(target_profile=output_profile,
+                                     sample=kwargs['sample'],
+                                     window_width=kwargs['w_width'],
+                                     window_height=kwargs['w_height'])
+
+            dsm = rio.open(kwargs['dsm_src'])
+            dtm = rio.open(kwargs['dtm_src'])
+            verify_matching_rasters(raster_a=dtm, raster_b=dsm)
+
+            for i, w in enumerate(windows, start=1):
+
+                print(f'processing window {i}/{len(windows)}', end="\r")
+                opg = get_opgaande_elementen(
+                    dsm=dsm,
+                    dtm=dtm,
+                    minimum_difference=kwargs['min_diff'],
+                    uniform_height=kwargs['u_height'],
+                    maaiveld_value=kwargs['mv_value'],
+                    window=w.rio_window,
+                )
+                dest.write(opg, 1,  window=w.rio_window)
+
+                u_values.update(np.unique(opg[opg != kwargs['mv_value']]))
+
+        if data_source == 'directory':
+
+            i = 0
+            for dsm_src, dtm_src in get_ahn_files(**kwargs):
+
+                dsm = rio.open(dsm_src)
+                dtm = rio.open(dtm_src)
+                verify_matching_rasters(raster_a=dtm, raster_b=dsm)
+
+                print(f'processing DSM-DTM file {i}/?', end="\r")
+                opg = get_opgaande_elementen(
+                    dsm=dsm,
+                    dtm=dtm,
+                    minimum_difference=kwargs['min_diff'],
+                    uniform_height=kwargs['u_height'],
+                    maaiveld_value=kwargs['mv_value'],
+                )
+                dest.write(opg, 1,  window=rio.windows.from_bounds(left=dsm.bounds.left,
+                                                       bottom=dsm.bounds.bottom,
+                                                       right=dsm.bounds.right,
+                                                       top=dsm.bounds.top,
+                                                       transform=output_profile['transform']))
+
+                u_values.update(np.unique(opg[opg != kwargs['mv_value']]))
+
+                i += 1
+
+        heights_table_to_file(os.path.join(kwargs['out_dir'],
                                        f'{kwargs["out_name"]}_heights.dbf'),
                           build_heights_table(list(u_values)))
 
@@ -95,7 +118,6 @@ if __name__ == '__main__':
 
     try:
         build_opgaande_elementenkaart(**vars(args))
-
         if args.filter_powerlines:
             filter_opgaande_elementen_under_powerlines(**vars(args))
     except (ValueError, AssertionError) as e:

sample/geslotenheid/ahn_opgaande_elemenden/ahn_read.py

@@ -160,10 +160,10 @@ def cap_dsm_to_dtm(dsm_array, dsm_nodata,
 def get_opgaande_elementen(
         dsm: rio.DatasetReader,
         dtm: rio.DatasetReader,
-        read_window: rio.windows.Window,
         minimum_difference: int = 3,
         maaiveld_value: int = 0,
         uniform_height: int = 0,
+        **kwargs
 ) -> np.array:
     """
     Read DSM and DTM for a certain bounding box. Return array with DSM-DTM differences
@@ -196,12 +196,12 @@ def get_opgaande_elementen(
     """
 
     # dtm_array = dtm.read(1)#, window=read_window)
-    dtm_array = fill_array(dtm.read(1, window=read_window), nodata=dtm.nodata)
-    dsm_array = dsm.read(1, window=read_window)
+    dtm_array = fill_array(dtm.read(1, **kwargs), nodata=dtm.nodata)
+    dsm_array = dsm.read(1, **kwargs)
 
     if np.all(dsm_array == dsm.nodata):
         # print('  All nodata, returning dummy array')
-        return np.multiply(np.ones(shape=(read_window.height, read_window.width)),
+        return np.multiply(np.ones(shape=dsm_array.shape),
                            maaiveld_value)
 
     verify = verify_dsm_exceeds_dtm(dsm_array=dsm_array, dsm_nodata=dsm.nodata,

sample/geslotenheid/ahn_opgaande_elemenden/output_profile.py

+import os
+import random
+from typing import Tuple
+
+import affine
+import numpy as np
+import pandas as pd
+import rasterio as rio
+from rasterio.coords import BoundingBox
+from rasterio.shutil import exists
+
+DEFAULT_BOUNDS = BoundingBox(left=0, bottom=300000, right=280000, top=625000)
+
+
+def build_output_profile(**kwargs) -> Tuple[dict, str]:
+
+    """
+    Build output profile
+
+    Parameters
+    ----------
+    kwargs
+
+    Returns
+    -------
+
+    """
+
+    if os.path.isfile(kwargs['dsm_src']):
+        dsm_sample_file = rio.open(kwargs['dsm_src'])
+        origin = 'file'
+    elif os.path.isdir(kwargs['dsm_src']):
+        dsm_sample_file = rio.open(os.path.join(kwargs['dsm_src'], os.listdir(kwargs['dsm_src'])[0]))
+        origin = 'directory'
+    else:
+        print('unknown source')
+
+    profile = rio.default_gtiff_profile.copy()
+    profile.update(dtype=np.uint16, nodata=0, crs=dsm_sample_file.crs)
+    profile.update(driver={'tif': "GTiff", 'flt': 'EHdr'}[kwargs['of']])
+
+    resolution = dsm_sample_file.res[0]
+
+    profile.update(transform=affine.Affine(a=resolution,
+                              b=0,
+                              c=dsm_sample_file.bounds.left if origin == "file" else DEFAULT_BOUNDS.left,
+                              d=0,
+                              e=resolution * -1,
+                              f=dsm_sample_file.bounds.top if origin == "file" else DEFAULT_BOUNDS.top
+    ))
+
+    profile.update(height=dsm_sample_file.shape[0] if origin == 'file' else int((DEFAULT_BOUNDS.top - DEFAULT_BOUNDS.bottom) / resolution),
+                   width=dsm_sample_file.shape[1] if origin == 'file' else int((DEFAULT_BOUNDS.right - DEFAULT_BOUNDS.left) / resolution),
+                   count=1,
+                   nodata=0)
+
+    return profile, origin
+
+
+def get_ahn_files(**kwargs):
+    """
+    retrieve list of DSM and DTM files. Match based on filename
+    Parameters
+    ----------
+    kwargs
+
+    Returns
+    -------
+
+    """
+
+
+    d = {}
+    for x in ['dsm', 'dtm']:
+        assert os.path.isdir(kwargs[f'{x}_src'])
+        d[x] = [os.path.join(kwargs[f'{x}_src'], f) for f in os.listdir(kwargs[f'{x}_src']) if exists(os.path.join(kwargs[f'{x}_src'], f))]
+    dsms = pd.Series(d['dsm']).sort_values()
+    dtms = pd.Series(d['dtm']).sort_values()
+    assert dsms.shape == dtms.shape
+
+    if kwargs['sample'] > 0:
+        return random.sample([(dsm, dtm) for dsm, dtm in zip(dsms, dtms)], kwargs['sample'])
+    else:
+        return [(dsm, dtm) for dsm, dtm in zip(dsms, dtms)]
+