Skip to content

Fix issue with Coord-sys. Set output coordsys directly from input dtb… #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
danpejobo wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Fix issue with Coord-sys. Set output coordsys directly from input dtb… #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
10d02a8
Fix issue with Coord-sys. Set output coordsys directly from input dtb…
Feb 9, 2024
b954d8b
Merge branch 'norwegian-geotechnical-institute:main' into main
danpejobo Mar 5, 2025
de3bb63
Skip points outside the raster bounds
Mar 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 3 additions & 11 deletions BegrensSkade.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
20 changes: 20 additions & 0 deletions BegrensSkadeLib.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand All @@ -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)
Expand All @@ -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

Expand Down