From 10d02a89f9fba7c7a8927b941a1dd2b3d41a9cae Mon Sep 17 00:00:00 2001 From: Daniel Date: Fri, 9 Feb 2024 22:52:52 +0100 Subject: [PATCH 1/2] Fix issue with Coord-sys. Set output coordsys directly from input dtb_raster. Already reprojected. --- BegrensSkade.py | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/BegrensSkade.py b/BegrensSkade.py index ccb90e2..a77724a 100644 --- a/BegrensSkade.py +++ b/BegrensSkade.py @@ -621,18 +621,10 @@ def mainBegrensSkade_ImpactMap( outBand.WriteArray(outData, 0, 0) outBand.SetNoDataValue(-99) - # georeference the image and set the projection + # Set the geotransform and projection from the input raster because the input + # is already in the correct projection outDs.SetGeoTransform(dataset.GetGeoTransform()) - #outDs.SetProjection(dataset.GetProjection()) - - projections = Utils.getProjections() - - output_def_desc = projections.get(str(output_proj)) - if not output_def_desc: - raise Exception(f"SRID: {output_proj} is not supported") - output_definition = output_def_desc["definition"]["data"] - - outDs.SetProjection(output_definition) + outDs.SetProjection(dataset.GetProjection()) outDs.FlushCache() del outData From de3bb6330c407450876d38639f99bbeb3864cfe1 Mon Sep 17 00:00:00 2001 From: Daniel Date: Thu, 6 Mar 2025 10:35:22 +0100 Subject: [PATCH 2/2] Skip points outside the raster bounds --- BegrensSkadeLib.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/BegrensSkadeLib.py b/BegrensSkadeLib.py index 56fcc18..73936b5 100644 --- a/BegrensSkadeLib.py +++ b/BegrensSkadeLib.py @@ -595,6 +595,15 @@ def get_buildings_with_dtb(features, dtb_filename, fieldNameFoundation=None, fie foundation = None structure = None status = None + + min_x = gt[0] + max_y = gt[3] + max_x = min_x + (gt[1] * src_ds.RasterXSize) + min_y = max_y + (gt[5] * src_ds.RasterYSize) # gt[5] is usually negative + + # Counter for points outside raster bounds + numPointsOutsideRaster = 0 + for buildingFeature in layer: @@ -615,6 +624,13 @@ def get_buildings_with_dtb(features, dtb_filename, fieldNameFoundation=None, fie for pnt in g.GetPoints(): if pnt: + x, y = pnt[0], pnt[1] + + # Skip points outside the raster bounds + if x < min_x or x > max_x or y < min_y or y > max_y: + numPointsOutsideRaster += 1 + continue + px = int((pnt[0] - gt[0]) / gt[1]) # x pixel py = int((pnt[1] - gt[3]) / gt[5]) # y pixel dtb = rb.ReadAsArray(px, py, 1, 1) @@ -640,6 +656,10 @@ def get_buildings_with_dtb(features, dtb_filename, fieldNameFoundation=None, fie Building(bid, corners, g.GetArea(), g.Length(), foundation=foundation, structure=structure, status=status, logger=logger)) bid += 1 + # Log the total number of skipped points **only once** + if numPointsOutsideRaster > 0: + logger.debug(f"Total points skipped due to falling outside raster bounds: {numPointsOutsideRaster}") + logger.debug(f"Number of skipped buildings: {numSkippedBuildings} of {numBuildings} - get_buildings_with_dtb") return input_buildings