From 886fc56d0dcf32b1a2c345e6403f21dc9dece0f2 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 12:47:39 -0400 Subject: [PATCH 01/21] Rewrote raster export (as-is case) --- app/raster_utils.py | 64 +++++++++++++++++++++++++-------------------- 1 file changed, 36 insertions(+), 28 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 91b48a25..217bc060 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -42,6 +42,9 @@ def export_raster(input, output, **opts): dem = asset_type in ['dsm', 'dtm'] + import time + now = time.time() + if crop_wkt is not None: with rasterio.open(input) as ds: crop = GEOSGeometry(crop_wkt) @@ -137,7 +140,10 @@ def export_raster(input, output, **opts): rescale = [md['statistics']['1']['min'], md['statistics']['1']['max']] ci = src.colorinterp - + alpha_index = None + if has_alpha_band(src): + alpha_index = src.colorinterp.index(ColorInterp.alpha) + 1 + if rgb and expression is None: # More than 4 bands? if len(ci) > 4: @@ -170,14 +176,29 @@ def export_raster(input, output, **opts): def process(arr, skip_rescale=False, skip_alpha=False, skip_type=False): + aidx = None + if alpha_index is not None: + aidx = indexes.index(alpha_index) + + # Alpha must (should) be the last band + if aidx != len(indexes) - 1: + # TODO remove + print("Alpha index is not the last one") + aidx = None + bands, w, h = arr.shape + if not skip_rescale and rescale is not None: - arr = linear_rescale(arr, in_range=rescale) - if not skip_alpha and not with_alpha: - arr[mask==0] = jpg_background + arr[:aidx,:,:] = linear_rescale(arr[:aidx,:,:], in_range=rescale) + if not skip_alpha and not with_alpha and aidx is not None: + background = arr[aidx]==0 + arr[:aidx, :, :][:, background] = jpg_background if not skip_type and rgb and arr.dtype != np.uint8: arr = arr.astype(np.uint8) - return arr + if not with_alpha and aidx is not None: + return arr[:aidx,:,:] + else: + return arr def update_rgb_colorinterp(dst): if with_alpha: @@ -217,9 +238,9 @@ def export_raster(input, output, **opts): height=height ) - def write_band(arr, dst, band_num): + def write_to(arr, dst): reproject(source=arr, - destination=rasterio.band(dst, band_num), + destination=dst, src_transform=win_transform, src_crs=src.crs, dst_transform=transform, @@ -229,8 +250,8 @@ def export_raster(input, output, **opts): else: # No reprojection needed - def write_band(arr, dst, band_num): - dst.write(arr, band_num) + def write_to(arr, dst): + dst.write(arr) if expression is not None: # Apply band math @@ -243,13 +264,8 @@ def export_raster(input, output, **opts): rgb_expr = expression.split(",") indexes = tuple([int(b.replace("b", "")) for b in bands_names]) - alpha_index = None - if has_alpha_band(src): - try: - alpha_index = src.colorinterp.index(ColorInterp.alpha) + 1 - indexes += (alpha_index, ) - except ValueError: - pass + if alpha_index is not None: + indexes += (alpha_index, ) data = reader.read(indexes=indexes, window=win, out_dtype=np.float32) arr = dict(zip(bands_names, data)) @@ -322,18 +338,8 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: - band_num = 1 - for idx in indexes: - ci = src.colorinterp[idx - 1] - arr = reader.read(idx, window=win) - - if ci == ColorInterp.alpha: - if with_alpha: - write_band(arr, dst, band_num) - band_num += 1 - else: - write_band(process(arr), dst, band_num) - band_num += 1 + arr = reader.read(indexes=indexes, window=win) + write_to(process(arr), dst) new_ci = [src.colorinterp[idx - 1] for idx in indexes] if not with_alpha: @@ -344,6 +350,8 @@ def export_raster(input, output, **opts): # Close warped vrt if vrt_options is not None: reader.close() + + logger.info(f"Finished in {time.time() - now}s") if kmz: subprocess.check_output(["gdal_translate", "-of", "KMLSUPEROVERLAY", From 04e9f269212c549b541baaaab7301c961f84f1a2 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 13:00:51 -0400 Subject: [PATCH 02/21] Fix reproject --- app/raster_utils.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 217bc060..6fd8c57d 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -162,11 +162,6 @@ def export_raster(input, output, **opts): indexes = (ci.index(ColorInterp.gray) + 1,) * 3 + \ (ci.index(ColorInterp.alpha) + 1, ) - if ColorInterp.alpha in ci: - mask = reader.read(ci.index(ColorInterp.alpha) + 1, window=win) - else: - mask = reader.dataset_mask(window=win) - cmap = None if color_map: try: @@ -240,7 +235,7 @@ def export_raster(input, output, **opts): def write_to(arr, dst): reproject(source=arr, - destination=dst, + destination=rasterio.band(dst, indexes), src_transform=win_transform, src_crs=src.crs, dst_transform=transform, From 8627da690b7b6b9c107010b63953cca57790d030 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 13:23:37 -0400 Subject: [PATCH 03/21] Windowed writes --- app/raster_utils.py | 47 ++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 42 insertions(+), 5 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 6fd8c57d..fb2c823e 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -29,6 +29,39 @@ def extension_for_export_format(export_format): } return extensions.get(export_format, export_format) +# Based on https://github.com/uav4geo/GeoDeep/blob/main/geodeep/slidingwindow.py +def compute_subwindows(window, max_window_size, overlap_pixels=0): + win_size_x = max_window_size + win_size_y = max_window_size + win_size_x = min(win_size_x, window.width) + win_size_y = min(win_size_y, window.height) + + step_size_x = win_size_x - overlap_pixels + step_size_y = win_size_y - overlap_pixels + + last_x = window.width - win_size_x + last_y = window.height - win_size_y + x_offsets = list(range(0, last_x + 1, step_size_x)) + y_offsets = list(range(0, last_y + 1, step_size_y)) + + if len(x_offsets) == 0 or x_offsets[-1] != last_x: + x_offsets.append(last_x) + if len(y_offsets) == 0 or y_offsets[-1] != last_y: + y_offsets.append(last_y) + + # Generate the list of windows + windows = [] + for x_offset in x_offsets: + for y_offset in y_offsets: + windows.append(Window( + x_offset, + y_offset, + win_size_x, + win_size_y, + )) + + return windows + def export_raster(input, output, **opts): epsg = opts.get('epsg') expression = opts.get('expression') @@ -233,7 +266,8 @@ def export_raster(input, output, **opts): height=height ) - def write_to(arr, dst): + def write_to(arr, dst, window=None): + # TODO reproject(source=arr, destination=rasterio.band(dst, indexes), src_transform=win_transform, @@ -245,8 +279,8 @@ def export_raster(input, output, **opts): else: # No reprojection needed - def write_to(arr, dst): - dst.write(arr) + def write_to(arr, dst, window=None): + dst.write(arr, window=window) if expression is not None: # Apply band math @@ -333,8 +367,11 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: - arr = reader.read(indexes=indexes, window=win) - write_to(process(arr), dst) + subwins = compute_subwindows(win, 100000) + for w in subwins: + logger.info(w) + arr = reader.read(indexes=indexes, window=w) + write_to(process(arr), dst, window=w) new_ci = [src.colorinterp[idx - 1] for idx in indexes] if not with_alpha: From 1b3b96cde0d33340b1dd11f46e4abd9025f30dc4 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 16:12:38 -0400 Subject: [PATCH 04/21] Fix windows calculation --- app/raster_utils.py | 56 ++++++++++++++++++++++++++++----------------- worker/tasks.py | 1 + 2 files changed, 36 insertions(+), 21 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index fb2c823e..434c0dc7 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -8,7 +8,7 @@ import numexpr as ne from django.contrib.gis.geos import GEOSGeometry from rasterio.enums import ColorInterp from rasterio.vrt import WarpedVRT -from rasterio.windows import Window, bounds as window_bounds +from rasterio.windows import Window, bounds as window_bounds, from_bounds from rio_tiler.utils import has_alpha_band, linear_rescale from rio_tiler.colormap import cmap as colormap, apply_cmap from rio_tiler.errors import InvalidColorMapName @@ -31,16 +31,19 @@ def extension_for_export_format(export_format): # Based on https://github.com/uav4geo/GeoDeep/blob/main/geodeep/slidingwindow.py def compute_subwindows(window, max_window_size, overlap_pixels=0): - win_size_x = max_window_size - win_size_y = max_window_size - win_size_x = min(win_size_x, window.width) - win_size_y = min(win_size_y, window.height) + col_off = int(window.col_off) + row_off = int(window.row_off) + width = int(window.width) + height = int(window.height) + + win_size_x = min(max_window_size, width) + win_size_y = min(max_window_size, height) step_size_x = win_size_x - overlap_pixels step_size_y = win_size_y - overlap_pixels - last_x = window.width - win_size_x - last_y = window.height - win_size_y + last_x = width - win_size_x - col_off + last_y = height - win_size_y - row_off x_offsets = list(range(0, last_x + 1, step_size_x)) y_offsets = list(range(0, last_y + 1, step_size_y)) @@ -54,8 +57,8 @@ def compute_subwindows(window, max_window_size, overlap_pixels=0): for x_offset in x_offsets: for y_offset in y_offsets: windows.append(Window( - x_offset, - y_offset, + x_offset + col_off, + y_offset + row_off, win_size_x, win_size_y, )) @@ -91,14 +94,14 @@ def export_raster(input, output, **opts): src = ds_src.dataset profile = src.meta.copy() win_bounds = None - win_transform = None + src_transform = None dst_width = src.width dst_height = src.height if vrt_options is None: reader = src win = Window(0, 0, dst_width, dst_height) - win_transform = src.transform + src_transform = src.transform else: reader = WarpedVRT(src, **vrt_options) dst_width = bounds[2] - bounds[0] @@ -106,7 +109,7 @@ def export_raster(input, output, **opts): win = Window(bounds[0], bounds[1], dst_width, dst_height) win_bounds = window_bounds(win, profile['transform']) - win_transform, width, height = calculate_default_transform( + src_transform, width, height = calculate_default_transform( src.crs, src.crs, src.width, src.height, left=win_bounds[0], bottom=win_bounds[1], @@ -116,7 +119,7 @@ def export_raster(input, output, **opts): dst_height=dst_height) profile.update( - transform=win_transform, + transform=src_transform, width=width, height=height ) @@ -241,13 +244,13 @@ def export_raster(input, output, **opts): if dem and rgb and profile.get('nodata') is not None: profile.update(nodata=None) - # Define write band function + # Define write to function # Reprojection needed? if src.crs is not None and epsg is not None and src.crs.to_epsg() != epsg: dst_crs = "EPSG:{}".format(epsg) if win_bounds is not None: - transform, width, height = calculate_default_transform( + dst_transform, width, height = calculate_default_transform( src.crs, dst_crs, src.width, src.height, left=win_bounds[0], bottom=win_bounds[1], @@ -256,12 +259,12 @@ def export_raster(input, output, **opts): dst_width=dst_width, dst_height=dst_height) else: - transform, width, height = calculate_default_transform( + dst_transform, width, height = calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) profile.update( crs=dst_crs, - transform=transform, + transform=dst_transform, width=width, height=height ) @@ -270,18 +273,29 @@ def export_raster(input, output, **opts): # TODO reproject(source=arr, destination=rasterio.band(dst, indexes), - src_transform=win_transform, + src_transform=src_transform, src_crs=src.crs, - dst_transform=transform, + dst_transform=dst_transform, dst_crs=dst_crs, resampling=Resampling.nearest, num_threads=4) else: # No reprojection needed + dst_transform = src_transform + def write_to(arr, dst, window=None): dst.write(arr, window=window) + def compute_dst_window(w): + dst_w = from_bounds( + *window_bounds(w, src_transform), + transform=dst_transform, + ) + + return Window(int(round(dst_w.col_off)), int(round(dst_w.row_off)), + int(round(dst_w.width)), int(round(dst_w.height))) + if expression is not None: # Apply band math if rgb: @@ -367,11 +381,11 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: - subwins = compute_subwindows(win, 100000) + subwins = compute_subwindows(win, 1024) for w in subwins: logger.info(w) arr = reader.read(indexes=indexes, window=w) - write_to(process(arr), dst, window=w) + write_to(process(arr), dst, window=compute_dst_window(w)) new_ci = [src.colorinterp[idx - 1] for idx in indexes] if not with_alpha: diff --git a/worker/tasks.py b/worker/tasks.py index 0b2fd353..02b5dd79 100644 --- a/worker/tasks.py +++ b/worker/tasks.py @@ -203,6 +203,7 @@ def export_raster(self, input, **opts): return result except Exception as e: + logger.error(traceback.format_exc()) logger.error(str(e)) return {'error': str(e)} From 40ccf94860945ce9a5e424b7c779efb96150d2e3 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 16:53:18 -0400 Subject: [PATCH 05/21] Fix compute subwindows --- app/raster_utils.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 434c0dc7..c2d71a73 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -42,10 +42,10 @@ def compute_subwindows(window, max_window_size, overlap_pixels=0): step_size_x = win_size_x - overlap_pixels step_size_y = win_size_y - overlap_pixels - last_x = width - win_size_x - col_off - last_y = height - win_size_y - row_off - x_offsets = list(range(0, last_x + 1, step_size_x)) - y_offsets = list(range(0, last_y + 1, step_size_y)) + last_x = col_off + width - win_size_x + last_y = row_off + height - win_size_y + x_offsets = list(range(col_off, last_x + 1, step_size_x)) + y_offsets = list(range(row_off, last_y + 1, step_size_y)) if len(x_offsets) == 0 or x_offsets[-1] != last_x: x_offsets.append(last_x) @@ -57,8 +57,8 @@ def compute_subwindows(window, max_window_size, overlap_pixels=0): for x_offset in x_offsets: for y_offset in y_offsets: windows.append(Window( - x_offset + col_off, - y_offset + row_off, + x_offset, + y_offset, win_size_x, win_size_y, )) @@ -293,7 +293,7 @@ def export_raster(input, output, **opts): transform=dst_transform, ) - return Window(int(round(dst_w.col_off)), int(round(dst_w.row_off)), + return Window(int(round(dst_w.col_off)) - win.col_off, int(round(dst_w.row_off)) - win.row_off, int(round(dst_w.width)), int(round(dst_w.height))) if expression is not None: @@ -381,12 +381,18 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: - subwins = compute_subwindows(win, 1024) + subwins = compute_subwindows(win, 2048) for w in subwins: logger.info(w) + logger.info(compute_dst_window(w)) + logger.info("====") arr = reader.read(indexes=indexes, window=w) write_to(process(arr), dst, window=compute_dst_window(w)) + + # arr = reader.read(indexes=indexes, window=win) + # write_to(process(arr), dst, window=compute_dst_window(win)) + new_ci = [src.colorinterp[idx - 1] for idx in indexes] if not with_alpha: new_ci = [ci for ci in new_ci if ci != ColorInterp.alpha] From fd63e3d32db760d13d434b9f5be21eebfaa42240 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 18:09:24 -0400 Subject: [PATCH 06/21] Improve compute_dst_window --- app/raster_utils.py | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index c2d71a73..8a2c7af5 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -16,7 +16,7 @@ from app.api.hsvblend import hsv_blend from app.api.hillshade import LightSource from app.geoutils import geom_transform_wkt_bbox from rio_tiler.io import COGReader -from rasterio.warp import calculate_default_transform, reproject, Resampling +from rasterio.warp import calculate_default_transform, reproject, Resampling, transform_bounds logger = logging.getLogger('app.logger') @@ -269,32 +269,46 @@ def export_raster(input, output, **opts): height=height ) - def write_to(arr, dst, window=None): - # TODO + def write_to(arr, dst, window): + # bnds = window_bounds(window, dst_transform) + + # _, width, height = calculate_default_transform( + # src.crs, src.crs, src.width, src.height, + # left=win_bounds[0], + # bottom=win_bounds[1], + # right=win_bounds[2], + # top=win_bounds[3], + # dst_width=dst_width, + # dst_height=dst_height) + dst_arr = np.zeros((arr.shape[0], int(window.height), int(window.width)), dtype=arr.dtype) reproject(source=arr, - destination=rasterio.band(dst, indexes), + destination=dst_arr, src_transform=src_transform, src_crs=src.crs, dst_transform=dst_transform, dst_crs=dst_crs, resampling=Resampling.nearest, num_threads=4) + + dst.write(dst_arr, window=window) else: # No reprojection needed dst_transform = src_transform + dst_crs = src.crs - def write_to(arr, dst, window=None): + def write_to(arr, dst, window): dst.write(arr, window=window) def compute_dst_window(w): + src_bounds = window_bounds(w, src_transform) + dst_bounds = transform_bounds(src.crs, dst_crs, *src_bounds, densify_pts=21) dst_w = from_bounds( - *window_bounds(w, src_transform), + *dst_bounds, transform=dst_transform, ) - - return Window(int(round(dst_w.col_off)) - win.col_off, int(round(dst_w.row_off)) - win.row_off, - int(round(dst_w.width)), int(round(dst_w.height))) + return Window(int(dst_w.col_off) - win.col_off, int(dst_w.row_off) - win.row_off, + int(dst_w.width), int(dst_w.height)) if expression is not None: # Apply band math @@ -381,6 +395,7 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: + logger.info(f"MAIN: {win}") subwins = compute_subwindows(win, 2048) for w in subwins: logger.info(w) From a17bddd0994f00bff1d701e86e45d218530c0724 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 18:46:25 -0400 Subject: [PATCH 07/21] No data initialization --- app/raster_utils.py | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 8a2c7af5..73cc5f2f 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -270,27 +270,19 @@ def export_raster(input, output, **opts): ) def write_to(arr, dst, window): - # bnds = window_bounds(window, dst_transform) - - # _, width, height = calculate_default_transform( - # src.crs, src.crs, src.width, src.height, - # left=win_bounds[0], - # bottom=win_bounds[1], - # right=win_bounds[2], - # top=win_bounds[3], - # dst_width=dst_width, - # dst_height=dst_height) - dst_arr = np.zeros((arr.shape[0], int(window.height), int(window.width)), dtype=arr.dtype) + # dst_arr = np.zeros((arr.shape[0], int(window.height), int(window.width)), dtype=arr.dtype) + # for band in range(arr.shape[0]): + current = dst.read(window=window) reproject(source=arr, - destination=dst_arr, + destination=current, src_transform=src_transform, src_crs=src.crs, dst_transform=dst_transform, dst_crs=dst_crs, resampling=Resampling.nearest, - num_threads=4) - - dst.write(dst_arr, window=window) + num_threads=4, + init_dest_nodata=False) + dst.write(current, window=window) else: # No reprojection needed @@ -395,8 +387,10 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: + pass + with rasterio.open(output_raster, 'r+', **profile) as dst: logger.info(f"MAIN: {win}") - subwins = compute_subwindows(win, 2048) + subwins = compute_subwindows(win, 512) for w in subwins: logger.info(w) logger.info(compute_dst_window(w)) From 0c981c5c7da2534d259e847634b49e281b20f66f Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 21:04:06 -0400 Subject: [PATCH 08/21] GDAL to the rescue --- app/raster_utils.py | 111 ++++++++++++++------------------------------ 1 file changed, 36 insertions(+), 75 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 73cc5f2f..a9d92eca 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -16,7 +16,7 @@ from app.api.hsvblend import hsv_blend from app.api.hillshade import LightSource from app.geoutils import geom_transform_wkt_bbox from rio_tiler.io import COGReader -from rasterio.warp import calculate_default_transform, reproject, Resampling, transform_bounds +from rasterio.warp import calculate_default_transform logger = logging.getLogger('app.logger') @@ -133,15 +133,29 @@ def export_raster(input, output, **opts): indexes = src.indexes output_raster = output jpg_background = 255 # white + reproject = src.crs is not None and epsg is not None and src.crs.to_epsg() != epsg + path_base, _ = os.path.splitext(output) # KMZ is special, we just export it as GeoTIFF # and then call GDAL to tile/package it kmz = export_format == "kmz" if kmz: export_format = "gtiff-rgb" - path_base, _ = os.path.splitext(output) output_raster = path_base + ".kmz.tif" + # JPG and PNG are exported to GeoTIFF only if reprojection is needed + jpg = export_format == "jpg" + png = export_format == "png" + if reproject: + if jpg: + export_format = 'gtiff-rgb' + path_base, _ = os.path.splitext(output) + output_raster = path_base + ".jpg.tif" + if png: + export_format = 'gtiff-rgb' + path_base, _ = os.path.splitext(output) + output_raster = path_base + ".png.tif" + if export_format == "jpg": driver = "JPEG" profile.update(quality=90) @@ -158,11 +172,18 @@ def export_raster(input, output, **opts): profile.update(BIGTIFF='IF_SAFER') band_count = 4 rgb = True + if jpg: + band_count = 3 + with_alpha = False else: compress = "DEFLATE" profile.update(BIGTIFF='IF_SAFER') band_count = src.count + if reproject: + path_base, _ = os.path.splitext(output_raster) + output_raster = path_base + ".base.tif" + if compress is not None: profile.update(compress=compress) profile.update(predictor=2 if compress == "DEFLATE" else 1) @@ -213,8 +234,6 @@ def export_raster(input, output, **opts): # Alpha must (should) be the last band if aidx != len(indexes) - 1: - # TODO remove - print("Alpha index is not the last one") aidx = None bands, w, h = arr.shape @@ -244,64 +263,6 @@ def export_raster(input, output, **opts): if dem and rgb and profile.get('nodata') is not None: profile.update(nodata=None) - # Define write to function - # Reprojection needed? - if src.crs is not None and epsg is not None and src.crs.to_epsg() != epsg: - dst_crs = "EPSG:{}".format(epsg) - - if win_bounds is not None: - dst_transform, width, height = calculate_default_transform( - src.crs, dst_crs, src.width, src.height, - left=win_bounds[0], - bottom=win_bounds[1], - right=win_bounds[2], - top=win_bounds[3], - dst_width=dst_width, - dst_height=dst_height) - else: - dst_transform, width, height = calculate_default_transform( - src.crs, dst_crs, src.width, src.height, *src.bounds) - - profile.update( - crs=dst_crs, - transform=dst_transform, - width=width, - height=height - ) - - def write_to(arr, dst, window): - # dst_arr = np.zeros((arr.shape[0], int(window.height), int(window.width)), dtype=arr.dtype) - # for band in range(arr.shape[0]): - current = dst.read(window=window) - reproject(source=arr, - destination=current, - src_transform=src_transform, - src_crs=src.crs, - dst_transform=dst_transform, - dst_crs=dst_crs, - resampling=Resampling.nearest, - num_threads=4, - init_dest_nodata=False) - dst.write(current, window=window) - - else: - # No reprojection needed - dst_transform = src_transform - dst_crs = src.crs - - def write_to(arr, dst, window): - dst.write(arr, window=window) - - def compute_dst_window(w): - src_bounds = window_bounds(w, src_transform) - dst_bounds = transform_bounds(src.crs, dst_crs, *src_bounds, densify_pts=21) - dst_w = from_bounds( - *dst_bounds, - transform=dst_transform, - ) - return Window(int(dst_w.col_off) - win.col_off, int(dst_w.row_off) - win.row_off, - int(dst_w.width), int(dst_w.height)) - if expression is not None: # Apply band math if rgb: @@ -387,20 +348,12 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: - pass - with rasterio.open(output_raster, 'r+', **profile) as dst: - logger.info(f"MAIN: {win}") subwins = compute_subwindows(win, 512) for w in subwins: - logger.info(w) - logger.info(compute_dst_window(w)) - logger.info("====") arr = reader.read(indexes=indexes, window=w) - write_to(process(arr), dst, window=compute_dst_window(w)) - - # arr = reader.read(indexes=indexes, window=win) - # write_to(process(arr), dst, window=compute_dst_window(win)) - + dst_w = Window(w.col_off - win.col_off, w.row_off - win.row_off, + w.width, w.height) + dst.write(process(arr), window=dst_w) new_ci = [src.colorinterp[idx - 1] for idx in indexes] if not with_alpha: @@ -412,9 +365,17 @@ def export_raster(input, output, **opts): if vrt_options is not None: reader.close() - logger.info(f"Finished in {time.time() - now}s") - if kmz: subprocess.check_output(["gdal_translate", "-of", "KMLSUPEROVERLAY", "-co", "Name={}".format(name), "-co", "FORMAT=AUTO", output_raster, output]) + elif reproject: + subprocess.check_output(["gdalwarp", "-r", "near", + "-multi", + "-wo", "NUM_THREADS=4", + "-t_srs", f"EPSG:{epsg}", + "--config", "GDAL_CACHEMAX", "25%", + output_raster, output]) + + logger.info(f"Finished in {time.time() - now}s") + From 65120015bad1405c13c85b7712ed5e34cbe4a787 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 21:09:13 -0400 Subject: [PATCH 09/21] Win size global --- app/raster_utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index a9d92eca..5c9ac088 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -128,6 +128,7 @@ def export_raster(input, output, **opts): driver = "GTiff" compress = None max_bands = 9999 + window_size = 512 with_alpha = True rgb = False indexes = src.indexes @@ -348,11 +349,10 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: - subwins = compute_subwindows(win, 512) + subwins = compute_subwindows(win, window_size) for w in subwins: arr = reader.read(indexes=indexes, window=w) - dst_w = Window(w.col_off - win.col_off, w.row_off - win.row_off, - w.width, w.height) + dst_w = Window(w.col_off - win.col_off, w.row_off - win.row_off, w.width, w.height) dst.write(process(arr), window=dst_w) new_ci = [src.colorinterp[idx - 1] for idx in indexes] From fe2740b398503aa6f46b9f496348c50bab1271df Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 16 Jul 2025 22:14:18 -0400 Subject: [PATCH 10/21] Started refactoring of process --- app/raster_utils.py | 110 +++++++++++++++++++++++--------------------- 1 file changed, 57 insertions(+), 53 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 5c9ac088..bbb5aaa5 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -56,12 +56,19 @@ def compute_subwindows(window, max_window_size, overlap_pixels=0): windows = [] for x_offset in x_offsets: for y_offset in y_offsets: - windows.append(Window( + w = Window( x_offset, y_offset, win_size_x, win_size_y, - )) + ) + dst_w = Window( + x_offset - window.col_off, + y_offset - window.row_off, + win_size_x, win_size_y + ) + + windows.append((w, dst_w)) return windows @@ -201,7 +208,9 @@ def export_raster(input, output, **opts): alpha_index = None if has_alpha_band(src): alpha_index = src.colorinterp.index(ColorInterp.alpha) + 1 - + + subwins = compute_subwindows(win, window_size) + if rgb and expression is None: # More than 4 bands? if len(ci) > 4: @@ -228,26 +237,23 @@ def export_raster(input, output, **opts): logger.warning("Invalid colormap {}".format(color_map)) - def process(arr, skip_rescale=False, skip_alpha=False, skip_type=False): - aidx = None - if alpha_index is not None: - aidx = indexes.index(alpha_index) - - # Alpha must (should) be the last band - if aidx != len(indexes) - 1: - aidx = None - bands, w, h = arr.shape + def process(arr, skip_rescale=False, skip_background=False, skip_type=False, mask=None): + # To include or exclude alpha + ax = -1 if alpha_index is not None else None if not skip_rescale and rescale is not None: - arr[:aidx,:,:] = linear_rescale(arr[:aidx,:,:], in_range=rescale) - if not skip_alpha and not with_alpha and aidx is not None: - background = arr[aidx]==0 - arr[:aidx, :, :][:, background] = jpg_background + arr[:ax,:,:] = linear_rescale(arr[:ax,:,:], in_range=rescale) + if not skip_background and not with_alpha and (ax is not None or mask is not None): + if mask is not None: + background = mask==0 + else: + background = arr[ax]==0 + arr[:ax, :, :][:, background] = jpg_background if not skip_type and rgb and arr.dtype != np.uint8: arr = arr.astype(np.uint8) - if not with_alpha and aidx is not None: - return arr[:aidx,:,:] + if not with_alpha and ax is not None: + return arr[:ax,:,:] else: return arr @@ -278,40 +284,40 @@ def export_raster(input, output, **opts): if alpha_index is not None: indexes += (alpha_index, ) - data = reader.read(indexes=indexes, window=win, out_dtype=np.float32) - arr = dict(zip(bands_names, data)) - arr = np.array([np.nan_to_num(ne.evaluate(bloc.strip(), local_dict=arr)) for bloc in rgb_expr]) - - # Set nodata values - index_band = arr[0] - if alpha_index is not None: - # -1 is the last band = alpha - index_band[data[-1] == 0] = -9999 - - # Remove infinity values - index_band[index_band>1e+30] = -9999 - index_band[index_band<-1e+30] = -9999 - - # Make sure this is float32 - arr = arr.astype(np.float32) - with rasterio.open(output_raster, 'w', **profile) as dst: - # Apply colormap? - if rgb and cmap is not None: - rgb_data, _ = apply_cmap(process(arr, skip_alpha=True), cmap) + for w, dst_w in subwins: + data = reader.read(indexes=indexes, window=w, out_dtype=np.float32) + arr = dict(zip(bands_names, data)) + arr = np.array([np.nan_to_num(ne.evaluate(bloc.strip(), local_dict=arr)) for bloc in rgb_expr]) - band_num = 1 - for b in rgb_data: - write_band(process(b, skip_rescale=True), dst, band_num) - band_num += 1 + # Set nodata values + index_band = arr[0] + if alpha_index is not None: + # -1 is the last band = alpha + index_band[data[-1] == 0] = -9999 - if with_alpha: - write_band(mask, dst, band_num) - - update_rgb_colorinterp(dst) - else: - # Raw - write_band(process(arr)[0], dst, 1) + # Remove infinity values + index_band[index_band>1e+30] = -9999 + index_band[index_band<-1e+30] = -9999 + + # Make sure this is float32 + arr = arr.astype(np.float32) + + # Apply colormap? + if rgb and cmap is not None: + rgb_data, mask = apply_cmap(process(arr, skip_background=True), cmap) + + logger.warning(rgb_data.shape) + dst.write(process(rgb_data, skip_rescale=True, mask=mask), window=dst_w, indexes=(1,2,3)) + + if with_alpha: + dst.write(mask, 4, window=dst_w) + + update_rgb_colorinterp(dst) + else: + # Raw + # write_band(process(arr)[0], dst, 1) + dst.write(process(arr)[0], window=dst_w) elif dem: # Apply hillshading, colormaps to elevation with rasterio.open(output_raster, 'w', **profile) as dst: @@ -328,7 +334,7 @@ def export_raster(input, output, **opts): # Apply colormap? if rgb and cmap is not None: - rgb_data, _ = apply_cmap(process(arr, skip_alpha=True), cmap) + rgb_data, _ = apply_cmap(process(arr, skip_background=True), cmap) arr = None if intensity is not None: @@ -349,10 +355,8 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: - subwins = compute_subwindows(win, window_size) - for w in subwins: + for w, dst_w in subwins: arr = reader.read(indexes=indexes, window=w) - dst_w = Window(w.col_off - win.col_off, w.row_off - win.row_off, w.width, w.height) dst.write(process(arr), window=dst_w) new_ci = [src.colorinterp[idx - 1] for idx in indexes] From 1bb588b5b95981578cb6f0d2028c68e3a8cf65ee Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 00:37:21 -0400 Subject: [PATCH 11/21] Refactored process() --- app/raster_utils.py | 38 ++++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index bbb5aaa5..88f87fee 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -237,23 +237,23 @@ def export_raster(input, output, **opts): logger.warning("Invalid colormap {}".format(color_map)) - def process(arr, skip_rescale=False, skip_background=False, skip_type=False, mask=None): - # To include or exclude alpha - ax = -1 if alpha_index is not None else None - + def process(arr, skip_rescale=False, skip_background=False, skip_type=False, mask=None, includes_alpha=True, drop_last_band=False): if not skip_rescale and rescale is not None: - arr[:ax,:,:] = linear_rescale(arr[:ax,:,:], in_range=rescale) - if not skip_background and not with_alpha and (ax is not None or mask is not None): + arr = linear_rescale(arr, in_range=rescale) + if not skip_background and (mask is not None or includes_alpha): if mask is not None: background = mask==0 + elif includes_alpha: + background = arr[-1]==0 + if includes_alpha: + arr[:-1, :, :][:, background] = jpg_background else: - background = arr[ax]==0 - arr[:ax, :, :][:, background] = jpg_background + arr[:, background] = jpg_background if not skip_type and rgb and arr.dtype != np.uint8: arr = arr.astype(np.uint8) - - if not with_alpha and ax is not None: - return arr[:ax,:,:] + + if drop_last_band: + return arr[:-1, :, :] else: return arr @@ -292,9 +292,11 @@ def export_raster(input, output, **opts): # Set nodata values index_band = arr[0] + mask = None if alpha_index is not None: # -1 is the last band = alpha - index_band[data[-1] == 0] = -9999 + mask = data[-1] != 0 + index_band[~mask] = -9999 # Remove infinity values index_band[index_band>1e+30] = -9999 @@ -305,13 +307,11 @@ def export_raster(input, output, **opts): # Apply colormap? if rgb and cmap is not None: - rgb_data, mask = apply_cmap(process(arr, skip_background=True), cmap) - - logger.warning(rgb_data.shape) - dst.write(process(rgb_data, skip_rescale=True, mask=mask), window=dst_w, indexes=(1,2,3)) + rgb_data, _ = apply_cmap(process(arr, skip_background=True, includes_alpha=False), cmap) + dst.write(process(rgb_data, skip_rescale=True, mask=mask, includes_alpha=False), window=dst_w, indexes=(1,2,3)) if with_alpha: - dst.write(mask, 4, window=dst_w) + dst.write(mask.astype(np.uint8) * 255, 4, window=dst_w) update_rgb_colorinterp(dst) else: @@ -357,7 +357,7 @@ def export_raster(input, output, **opts): with rasterio.open(output_raster, 'w', **profile) as dst: for w, dst_w in subwins: arr = reader.read(indexes=indexes, window=w) - dst.write(process(arr), window=dst_w) + dst.write(process(arr, drop_last_band=not with_alpha), window=dst_w) new_ci = [src.colorinterp[idx - 1] for idx in indexes] if not with_alpha: @@ -374,6 +374,8 @@ def export_raster(input, output, **opts): "-co", "Name={}".format(name), "-co", "FORMAT=AUTO", output_raster, output]) elif reproject: + # TODO: we could use A VRT for better compression + # https://gdal.org/en/stable/programs/gdalwarp.html#compressed-output subprocess.check_output(["gdalwarp", "-r", "near", "-multi", "-wo", "NUM_THREADS=4", From 3e7351ae31ad4b601072f054b5294f1f42122937 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 10:16:35 -0400 Subject: [PATCH 12/21] Fix compression --- app/raster_utils.py | 46 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 37 insertions(+), 9 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 88f87fee..d5303e36 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -138,6 +138,7 @@ def export_raster(input, output, **opts): window_size = 512 with_alpha = True rgb = False + bigtiff = False indexes = src.indexes output_raster = output jpg_background = 255 # white @@ -176,23 +177,26 @@ def export_raster(input, output, **opts): rgb = True elif export_format == "gtiff-rgb": compress = "JPEG" + bigtiff = True profile.update(jpeg_quality=90) - profile.update(BIGTIFF='IF_SAFER') band_count = 4 rgb = True if jpg: band_count = 3 with_alpha = False else: + bigtiff = True compress = "DEFLATE" - profile.update(BIGTIFF='IF_SAFER') band_count = src.count + if bigtiff: + profile.update(BIGTIFF='IF_SAFER') + if reproject: path_base, _ = os.path.splitext(output_raster) output_raster = path_base + ".base.tif" - if compress is not None: + if compress is not None and not reproject: profile.update(compress=compress) profile.update(predictor=2 if compress == "DEFLATE" else 1) @@ -374,14 +378,38 @@ def export_raster(input, output, **opts): "-co", "Name={}".format(name), "-co", "FORMAT=AUTO", output_raster, output]) elif reproject: - # TODO: we could use A VRT for better compression - # https://gdal.org/en/stable/programs/gdalwarp.html#compressed-output + output_vrt = path_base + ".vrt" + subprocess.check_output(["gdalwarp", "-r", "near", - "-multi", - "-wo", "NUM_THREADS=4", + "-of", "VRT", "-t_srs", f"EPSG:{epsg}", - "--config", "GDAL_CACHEMAX", "25%", - output_raster, output]) + output_raster, output_vrt]) + # subprocess.check_output(["gdalwarp", "-r", "near", + # "-multi", + # "-wo", "NUM_THREADS=4", + # "-t_srs", f"EPSG:{epsg}", + # "--config", "GDAL_CACHEMAX", "25%", + # output_raster, output]) + gt_args = ["-r", "nearest", "--config", "GDAL_CACHEMAX", "25%"] + if bigtiff and not jpg and not png: + gt_args += ["-co", "BIGTIFF=IF_SAFER", + "-co", "BLOCKXSIZE=512", + "-co", "BLOCKYSIZE=512", + "-co", "NUM_THREADS=4",] + + if compress and not png: + if jpg: + gt_args += ["-co", "QUALITY=90"] + else: + gt_args += ["-co", f"COMPRESS={compress}", + "-co", "PREDICTOR=2"] + + subprocess.check_output(["gdal_translate"] + + gt_args + + [output_vrt, output]) + + if os.path.isfile(output_raster): + os.unlink(output_raster) logger.info(f"Finished in {time.time() - now}s") From 881960e384f94c43da475d833babfe48f09f5d93 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 10:47:25 -0400 Subject: [PATCH 13/21] Raster DEM exports --- app/raster_utils.py | 70 +++++++++++++++++++++------------------------ 1 file changed, 32 insertions(+), 38 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index d5303e36..5f961f6b 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -134,7 +134,6 @@ def export_raster(input, output, **opts): # Output format driver = "GTiff" compress = None - max_bands = 9999 window_size = 512 with_alpha = True rgb = False @@ -320,42 +319,43 @@ def export_raster(input, output, **opts): update_rgb_colorinterp(dst) else: # Raw - # write_band(process(arr)[0], dst, 1) - dst.write(process(arr)[0], window=dst_w) + dst.write(process(arr), window=dst_w) elif dem: # Apply hillshading, colormaps to elevation with rasterio.open(output_raster, 'w', **profile) as dst: - arr = reader.read(window=win) - - intensity = None - if hillshade is not None and hillshade > 0: - delta_scale = ZOOM_EXTRA_LEVELS ** 2 - dx = src.meta["transform"][0] * delta_scale - dy = src.meta["transform"][4] * delta_scale - ls = LightSource(azdeg=315, altdeg=45) - intensity = ls.hillshade(arr[0], dx=dx, dy=dy, vert_exag=hillshade) - intensity = intensity * 255.0 - - # Apply colormap? - if rgb and cmap is not None: - rgb_data, _ = apply_cmap(process(arr, skip_background=True), cmap) - arr = None - - if intensity is not None: - rgb_data = hsv_blend(rgb_data, intensity) - - band_num = 1 - for b in rgb_data: - write_band(process(b, skip_rescale=True), dst, band_num) - band_num += 1 + for w, dst_w in subwins: + arr = reader.read(window=w)[:1] - if with_alpha: - write_band(mask, dst, band_num) + # Apply colormap? + if rgb and cmap is not None: + nodata = profile.get('nodata') + if nodata is None: + nodata = -9999 + + mask = arr[0] != nodata - update_rgb_colorinterp(dst) - else: - # Raw - write_band(process(arr)[0], dst, 1) + intensity = None + if hillshade is not None and hillshade > 0: + delta_scale = ZOOM_EXTRA_LEVELS ** 2 + dx = src.meta["transform"][0] * delta_scale + dy = src.meta["transform"][4] * delta_scale + ls = LightSource(azdeg=315, altdeg=45) + intensity = ls.hillshade(arr[0], dx=dx, dy=dy, vert_exag=hillshade) + intensity = intensity * 255.0 + + rgb_data, _ = apply_cmap(process(arr, skip_background=True, includes_alpha=False), cmap) + + if intensity is not None: + rgb_data = hsv_blend(rgb_data, intensity) + + dst.write(process(rgb_data, skip_rescale=True, mask=mask, includes_alpha=False), window=dst_w, indexes=(1,2,3)) + if with_alpha: + dst.write(mask.astype(np.uint8) * 255, 4, window=dst_w) + + update_rgb_colorinterp(dst) + else: + # Raw + dst.write(process(arr), window=dst_w) else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: @@ -384,12 +384,6 @@ def export_raster(input, output, **opts): "-of", "VRT", "-t_srs", f"EPSG:{epsg}", output_raster, output_vrt]) - # subprocess.check_output(["gdalwarp", "-r", "near", - # "-multi", - # "-wo", "NUM_THREADS=4", - # "-t_srs", f"EPSG:{epsg}", - # "--config", "GDAL_CACHEMAX", "25%", - # output_raster, output]) gt_args = ["-r", "nearest", "--config", "GDAL_CACHEMAX", "25%"] if bigtiff and not jpg and not png: gt_args += ["-co", "BIGTIFF=IF_SAFER", From 1518a80060760799b49177ef7b42de82d5efb621 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 12:57:27 -0400 Subject: [PATCH 14/21] Fix padding --- app/raster_utils.py | 31 ++++++++++++++++++++++++------- 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 5f961f6b..5afb9c36 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -72,6 +72,9 @@ def compute_subwindows(window, max_window_size, overlap_pixels=0): return windows +def padded_window(w, pad): + return Window(w.col_off - pad, w.row_off - pad, w.width + pad * 2, w.height + pad * 2) + def export_raster(input, output, **opts): epsg = opts.get('epsg') expression = opts.get('expression') @@ -324,15 +327,25 @@ def export_raster(input, output, **opts): # Apply hillshading, colormaps to elevation with rasterio.open(output_raster, 'w', **profile) as dst: for w, dst_w in subwins: - arr = reader.read(window=w)[:1] - # Apply colormap? if rgb and cmap is not None: nodata = profile.get('nodata') if nodata is None: nodata = -9999 - - mask = arr[0] != nodata + + pad = 16 + elevation = reader.read(window=padded_window(w, pad), boundless=True, fill_value=nodata, out_shape=( + 1, + window_size + pad * 2, + window_size + pad * 2, + ), resampling=rasterio.enums.Resampling.bilinear)[:1][0] + + elevation[0:pad, 0:pad] = nodata + elevation[pad+window_size:pad*2+window_size, 0:pad] = nodata + elevation[0:pad, pad+window_size:pad*2+window_size] = nodata + elevation[pad+window_size:pad*2+window_size, pad+window_size:pad*2+window_size] = nodata + + mask = elevation != nodata intensity = None if hillshade is not None and hillshade > 0: @@ -340,14 +353,17 @@ def export_raster(input, output, **opts): dx = src.meta["transform"][0] * delta_scale dy = src.meta["transform"][4] * delta_scale ls = LightSource(azdeg=315, altdeg=45) - intensity = ls.hillshade(arr[0], dx=dx, dy=dy, vert_exag=hillshade) + + intensity = ls.hillshade(elevation, dx=dx, dy=dy, vert_exag=hillshade) + intensity = intensity[pad:pad+window_size, pad:pad+window_size] intensity = intensity * 255.0 - rgb_data, _ = apply_cmap(process(arr, skip_background=True, includes_alpha=False), cmap) + rgb_data, _ = apply_cmap(process(elevation[pad:window_size+pad, pad:window_size+pad][np.newaxis,:], skip_background=True, includes_alpha=False), cmap) if intensity is not None: rgb_data = hsv_blend(rgb_data, intensity) + mask = mask[pad:window_size+pad, pad:window_size+pad] dst.write(process(rgb_data, skip_rescale=True, mask=mask, includes_alpha=False), window=dst_w, indexes=(1,2,3)) if with_alpha: dst.write(mask.astype(np.uint8) * 255, 4, window=dst_w) @@ -355,7 +371,8 @@ def export_raster(input, output, **opts): update_rgb_colorinterp(dst) else: # Raw - dst.write(process(arr), window=dst_w) + arr = reader.read(window=w)[:1] + dst.write(process(arr), window=dst_w) else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: From 0deef8e65aff8d4011e352d62eb9457c13720886 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 13:33:00 -0400 Subject: [PATCH 15/21] Fix multispectral rgb export --- app/raster_utils.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 5afb9c36..64847d23 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -245,7 +245,10 @@ def export_raster(input, output, **opts): def process(arr, skip_rescale=False, skip_background=False, skip_type=False, mask=None, includes_alpha=True, drop_last_band=False): if not skip_rescale and rescale is not None: - arr = linear_rescale(arr, in_range=rescale) + if includes_alpha: + arr[:-1, :, :] = linear_rescale(arr[:-1, :, :], in_range=rescale) + else: + arr = linear_rescale(arr, in_range=rescale) if not skip_background and (mask is not None or includes_alpha): if mask is not None: background = mask==0 @@ -256,6 +259,9 @@ def export_raster(input, output, **opts): else: arr[:, background] = jpg_background if not skip_type and rgb and arr.dtype != np.uint8: + if includes_alpha: + arr[-1][arr[-1] > 255] = 255 + arr = arr.astype(np.uint8) if drop_last_band: From 617aac6a929881b53da79db9d8e64073ad5ae007 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 13:34:37 -0400 Subject: [PATCH 16/21] Bump version --- package.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package.json b/package.json index b14ce21b..e58b0a2e 100644 --- a/package.json +++ b/package.json @@ -1,6 +1,6 @@ { "name": "WebODM", - "version": "2.8.4", + "version": "2.9.0", "description": "User-friendly, extendable application and API for processing aerial imagery.", "main": "index.js", "scripts": { From ba2e23cca21ebff35be4edb422df6f120cbe0b78 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 14:08:09 -0400 Subject: [PATCH 17/21] Add export progress --- app/raster_utils.py | 35 +++++++++++++++---- .../app/js/components/ExportAssetPanel.jsx | 11 +++--- worker/tasks.py | 7 ++-- 3 files changed, 40 insertions(+), 13 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 64847d23..bd56bda2 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -5,6 +5,7 @@ import os import subprocess import numpy as np import numexpr as ne +import time from django.contrib.gis.geos import GEOSGeometry from rasterio.enums import ColorInterp from rasterio.vrt import WarpedVRT @@ -75,7 +76,14 @@ def compute_subwindows(window, max_window_size, overlap_pixels=0): def padded_window(w, pad): return Window(w.col_off - pad, w.row_off - pad, w.width + pad * 2, w.height + pad * 2) -def export_raster(input, output, **opts): +def export_raster(input, output, progress_callback=None, **opts): + current_progress = 0 + def p(text, perc=0): + nonlocal current_progress + current_progress += perc + if progress_callback is not None: + progress_callback(text, current_progress) + epsg = opts.get('epsg') expression = opts.get('expression') export_format = opts.get('format') @@ -87,8 +95,6 @@ def export_raster(input, output, **opts): crop_wkt = opts.get('crop') dem = asset_type in ['dsm', 'dtm'] - - import time now = time.time() if crop_wkt is not None: @@ -282,6 +288,11 @@ def export_raster(input, output, **opts): if dem and rgb and profile.get('nodata') is not None: profile.update(nodata=None) + + post_perc = 20 if reproject or kmz else 0 + num_wins = len(subwins) + progress_per_win = (100 - post_perc) / num_wins if num_wins > 0 else 0 + if expression is not None: # Apply band math if rgb: @@ -297,7 +308,9 @@ def export_raster(input, output, **opts): indexes += (alpha_index, ) with rasterio.open(output_raster, 'w', **profile) as dst: - for w, dst_w in subwins: + for idx, (w, dst_w) in enumerate(subwins): + p(f"Processing tile {idx}/{num_wins}", progress_per_win) + data = reader.read(indexes=indexes, window=w, out_dtype=np.float32) arr = dict(zip(bands_names, data)) arr = np.array([np.nan_to_num(ne.evaluate(bloc.strip(), local_dict=arr)) for bloc in rgb_expr]) @@ -332,7 +345,9 @@ def export_raster(input, output, **opts): elif dem: # Apply hillshading, colormaps to elevation with rasterio.open(output_raster, 'w', **profile) as dst: - for w, dst_w in subwins: + for idx, (w, dst_w) in enumerate(subwins): + p(f"Processing tile {idx}/{num_wins}", progress_per_win) + # Apply colormap? if rgb and cmap is not None: nodata = profile.get('nodata') @@ -382,7 +397,9 @@ def export_raster(input, output, **opts): else: # Copy bands as-is with rasterio.open(output_raster, 'w', **profile) as dst: - for w, dst_w in subwins: + for idx, (w, dst_w) in enumerate(subwins): + p(f"Processing tile {idx}/{num_wins}", progress_per_win) + arr = reader.read(indexes=indexes, window=w) dst.write(process(arr, drop_last_band=not with_alpha), window=dst_w) @@ -400,6 +417,8 @@ def export_raster(input, output, **opts): subprocess.check_output(["gdal_translate", "-of", "KMLSUPEROVERLAY", "-co", "Name={}".format(name), "-co", "FORMAT=AUTO", output_raster, output]) + p("Finalizing", post_perc) + elif reproject: output_vrt = path_base + ".vrt" @@ -428,5 +447,7 @@ def export_raster(input, output, **opts): if os.path.isfile(output_raster): os.unlink(output_raster) - logger.info(f"Finished in {time.time() - now}s") + p("Finalizing", post_perc) + + logger.info(f"Exported {output} in {round(time.time() - now, 2)}s") diff --git a/app/static/app/js/components/ExportAssetPanel.jsx b/app/static/app/js/components/ExportAssetPanel.jsx index 2638b8ed..84424f9d 100644 --- a/app/static/app/js/components/ExportAssetPanel.jsx +++ b/app/static/app/js/components/ExportAssetPanel.jsx @@ -76,7 +76,8 @@ export default class ExportAssetPanel extends React.Component { epsg: this.props.task.epsg || null, customEpsg: Storage.getItem("last_export_custom_epsg") || "4326", resample: 0, - exporting: false + exporting: false, + progress: null } } @@ -132,7 +133,7 @@ export default class ExportAssetPanel extends React.Component { if (typeof cb !== 'function') cb = undefined; const { task } = this.props; - this.setState({exporting: true, error: ""}); + this.setState({exporting: true, error: "", progress: null}); const data = this.getExportParams(format); if (this.state.epsg === "custom") Storage.setItem("last_export_custom_epsg", data.epsg); @@ -152,6 +153,8 @@ export default class ExportAssetPanel extends React.Component { Workers.downloadFile(result.celery_task_id, result.filename); if (cb !== undefined) cb(); } + }, (_, progress) => { + this.setState({progress}); }); }else if (result.url){ // Simple download @@ -182,7 +185,7 @@ export default class ExportAssetPanel extends React.Component { } render(){ - const {epsg, customEpsg, exporting, format, resample } = this.state; + const {epsg, customEpsg, exporting, format, resample, progress } = this.state; const { exportFormats } = this.props; const utmEPSG = this.props.task.epsg; @@ -233,7 +236,7 @@ export default class ExportAssetPanel extends React.Component {
    diff --git a/worker/tasks.py b/worker/tasks.py index 02b5dd79..1cf5d8fc 100644 --- a/worker/tasks.py +++ b/worker/tasks.py @@ -195,7 +195,10 @@ def export_raster(self, input, **opts): try: logger.info("Exporting raster {} with options: {}".format(input, json.dumps(opts))) tmpfile = tempfile.mktemp('_raster.{}'.format(extension_for_export_format(opts.get('format', 'gtiff'))), dir=settings.MEDIA_TMP) - export_raster_sync(input, tmpfile, **opts) + def progress_callback(status, perc): + self.update_state(state="PROGRESS", meta={"status": status, "progress": perc}) + + export_raster_sync(input, tmpfile, progress_callback=progress_callback, **opts) result = {'file': tmpfile} if settings.TESTING: @@ -203,7 +206,7 @@ def export_raster(self, input, **opts): return result except Exception as e: - logger.error(traceback.format_exc()) + # logger.error(traceback.format_exc()) logger.error(str(e)) return {'error': str(e)} From 4d44338a50ca7178b06c739a6554acceacb526b7 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 15:07:11 -0400 Subject: [PATCH 18/21] Boundless reads with VRT via GDAL --- app/raster_utils.py | 84 +++++++++++++++++---------------------------- worker/tasks.py | 2 +- 2 files changed, 32 insertions(+), 54 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index bd56bda2..96d53d0b 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -8,16 +8,14 @@ import numexpr as ne import time from django.contrib.gis.geos import GEOSGeometry from rasterio.enums import ColorInterp -from rasterio.vrt import WarpedVRT from rasterio.windows import Window, bounds as window_bounds, from_bounds from rio_tiler.utils import has_alpha_band, linear_rescale from rio_tiler.colormap import cmap as colormap, apply_cmap from rio_tiler.errors import InvalidColorMapName from app.api.hsvblend import hsv_blend from app.api.hillshade import LightSource -from app.geoutils import geom_transform_wkt_bbox -from rio_tiler.io import COGReader from rasterio.warp import calculate_default_transform +from rio_tiler.io import COGReader logger = logging.getLogger('app.logger') @@ -77,6 +75,8 @@ def padded_window(w, pad): return Window(w.col_off - pad, w.row_off - pad, w.width + pad * 2, w.height + pad * 2) def export_raster(input, output, progress_callback=None, **opts): + now = time.time() + current_progress = 0 def p(text, perc=0): nonlocal current_progress @@ -95,50 +95,33 @@ def export_raster(input, output, progress_callback=None, **opts): crop_wkt = opts.get('crop') dem = asset_type in ['dsm', 'dtm'] - now = time.time() + path_base, _ = os.path.splitext(output) + resampling = 'nearest' + if dem: + resampling = 'bilinear' if crop_wkt is not None: - with rasterio.open(input) as ds: - crop = GEOSGeometry(crop_wkt) - crop.srid = 4326 - cutline, bounds = geom_transform_wkt_bbox(crop, ds, 'raster') - vrt_options = {'cutline': cutline, 'nodata': 0} - else: - vrt_options = None + crop = GEOSGeometry(crop_wkt) + crop.srid = 4326 + + crop_geojson = os.path.join(path_base, "crop.geojson") + raster_vrt = os.path.join(path_base, "raster.vrt") + + os.makedirs(os.path.dirname(crop_geojson), exist_ok=True) + with open(crop_geojson, "w", encoding="utf-8") as f: + f.write(crop.geojson) + + subprocess.check_output(["gdalwarp", "-cutline", crop_geojson, + '--config', 'GDALWARP_DENSIFY_CUTLINE', 'NO', + '-crop_to_cutline', '-of', 'VRT', '-r', resampling, + input, raster_vrt]) + + input = raster_vrt - with COGReader(input, vrt_options=vrt_options) as ds_src: + with COGReader(input) as ds_src: src = ds_src.dataset profile = src.meta.copy() - win_bounds = None - src_transform = None - dst_width = src.width - dst_height = src.height - - if vrt_options is None: - reader = src - win = Window(0, 0, dst_width, dst_height) - src_transform = src.transform - else: - reader = WarpedVRT(src, **vrt_options) - dst_width = bounds[2] - bounds[0] - dst_height = bounds[3] - bounds[1] - win = Window(bounds[0], bounds[1], dst_width, dst_height) - win_bounds = window_bounds(win, profile['transform']) - - src_transform, width, height = calculate_default_transform( - src.crs, src.crs, src.width, src.height, - left=win_bounds[0], - bottom=win_bounds[1], - right=win_bounds[2], - top=win_bounds[3], - dst_width=dst_width, - dst_height=dst_height) - - profile.update( - transform=src_transform, - width=width, - height=height - ) + win = Window(0, 0, src.width, src.height) # Output format driver = "GTiff" @@ -151,7 +134,6 @@ def export_raster(input, output, progress_callback=None, **opts): output_raster = output jpg_background = 255 # white reproject = src.crs is not None and epsg is not None and src.crs.to_epsg() != epsg - path_base, _ = os.path.splitext(output) # KMZ is special, we just export it as GeoTIFF # and then call GDAL to tile/package it @@ -213,7 +195,7 @@ def export_raster(input, output, progress_callback=None, **opts): nodata = None if asset_type == 'orthophoto': nodata = 0 - md = ds_src.metadata(pmin=2.0, pmax=98.0, hist_options={"bins": 255}, nodata=nodata, vrt_options=vrt_options) + md = ds_src.metadata(pmin=2.0, pmax=98.0, hist_options={"bins": 255}, nodata=nodata) rescale = [md['statistics']['1']['min'], md['statistics']['1']['max']] ci = src.colorinterp @@ -311,7 +293,7 @@ def export_raster(input, output, progress_callback=None, **opts): for idx, (w, dst_w) in enumerate(subwins): p(f"Processing tile {idx}/{num_wins}", progress_per_win) - data = reader.read(indexes=indexes, window=w, out_dtype=np.float32) + data = src.read(indexes=indexes, window=w, out_dtype=np.float32) arr = dict(zip(bands_names, data)) arr = np.array([np.nan_to_num(ne.evaluate(bloc.strip(), local_dict=arr)) for bloc in rgb_expr]) @@ -355,7 +337,7 @@ def export_raster(input, output, progress_callback=None, **opts): nodata = -9999 pad = 16 - elevation = reader.read(window=padded_window(w, pad), boundless=True, fill_value=nodata, out_shape=( + elevation = src.read(window=padded_window(w, pad), boundless=True, fill_value=nodata, out_shape=( 1, window_size + pad * 2, window_size + pad * 2, @@ -392,7 +374,7 @@ def export_raster(input, output, progress_callback=None, **opts): update_rgb_colorinterp(dst) else: # Raw - arr = reader.read(window=w)[:1] + arr = src.read(window=w)[:1] dst.write(process(arr), window=dst_w) else: # Copy bands as-is @@ -400,7 +382,7 @@ def export_raster(input, output, progress_callback=None, **opts): for idx, (w, dst_w) in enumerate(subwins): p(f"Processing tile {idx}/{num_wins}", progress_per_win) - arr = reader.read(indexes=indexes, window=w) + arr = src.read(indexes=indexes, window=w) dst.write(process(arr, drop_last_band=not with_alpha), window=dst_w) new_ci = [src.colorinterp[idx - 1] for idx in indexes] @@ -409,10 +391,6 @@ def export_raster(input, output, progress_callback=None, **opts): dst.colorinterp = new_ci - # Close warped vrt - if vrt_options is not None: - reader.close() - if kmz: subprocess.check_output(["gdal_translate", "-of", "KMLSUPEROVERLAY", "-co", "Name={}".format(name), @@ -426,7 +404,7 @@ def export_raster(input, output, progress_callback=None, **opts): "-of", "VRT", "-t_srs", f"EPSG:{epsg}", output_raster, output_vrt]) - gt_args = ["-r", "nearest", "--config", "GDAL_CACHEMAX", "25%"] + gt_args = ["-r", resampling, "--config", "GDAL_CACHEMAX", "25%"] if bigtiff and not jpg and not png: gt_args += ["-co", "BIGTIFF=IF_SAFER", "-co", "BLOCKXSIZE=512", diff --git a/worker/tasks.py b/worker/tasks.py index 1cf5d8fc..fac1a18d 100644 --- a/worker/tasks.py +++ b/worker/tasks.py @@ -206,7 +206,7 @@ def export_raster(self, input, **opts): return result except Exception as e: - # logger.error(traceback.format_exc()) + logger.error(traceback.format_exc()) # TODO uncomment logger.error(str(e)) return {'error': str(e)} From f6d7a3a3aed409993f41f39c9e8c0cc5fa3eae1e Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 15:11:34 -0400 Subject: [PATCH 19/21] Cleanup --- worker/tasks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/worker/tasks.py b/worker/tasks.py index fac1a18d..1cf5d8fc 100644 --- a/worker/tasks.py +++ b/worker/tasks.py @@ -206,7 +206,7 @@ def export_raster(self, input, **opts): return result except Exception as e: - logger.error(traceback.format_exc()) # TODO uncomment + # logger.error(traceback.format_exc()) logger.error(str(e)) return {'error': str(e)} From 9d1c664b28fd6e3057c5e5f4473ba572f8ffdcdc Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 15:14:16 -0400 Subject: [PATCH 20/21] Cleanup imports --- app/raster_utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 96d53d0b..891c0264 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -8,13 +8,12 @@ import numexpr as ne import time from django.contrib.gis.geos import GEOSGeometry from rasterio.enums import ColorInterp -from rasterio.windows import Window, bounds as window_bounds, from_bounds +from rasterio.windows import Window from rio_tiler.utils import has_alpha_band, linear_rescale from rio_tiler.colormap import cmap as colormap, apply_cmap from rio_tiler.errors import InvalidColorMapName from app.api.hsvblend import hsv_blend from app.api.hillshade import LightSource -from rasterio.warp import calculate_default_transform from rio_tiler.io import COGReader logger = logging.getLogger('app.logger') From bb13ca32ee06c2e7a09150ab45baf3e518b5511d Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 17 Jul 2025 15:19:53 -0400 Subject: [PATCH 21/21] Pass interpolation to gdalwarp --- app/raster_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/app/raster_utils.py b/app/raster_utils.py index 891c0264..e3814ca2 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -399,7 +399,7 @@ def export_raster(input, output, progress_callback=None, **opts): elif reproject: output_vrt = path_base + ".vrt" - subprocess.check_output(["gdalwarp", "-r", "near", + subprocess.check_output(["gdalwarp", "-r", "near" if resampling == "nearest" else resampling, "-of", "VRT", "-t_srs", f"EPSG:{epsg}", output_raster, output_vrt])