diff --git a/BegrensSkade.py b/BegrensSkade.py index 97c17de..0b59be8 100644 --- a/BegrensSkade.py +++ b/BegrensSkade.py @@ -603,18 +603,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 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