Skip to content

[#300] Fix for raster initialization code used in sinmap python scripts #301

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

Merged
dtarb merged 3 commits into from
May 19, 2026
Merged

[#300] Fix for raster initialization code used in sinmap python scripts #301

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8901198
fix: Fix numpy data type attribute error
pkdash May 18, 2026
51e6c74
fix: Improve raster projection handling in initialize_output_raster_f…
pkdash May 18, 2026
4b94b84
docs: add comments to clarify raster projection setting code logic
pkdash May 19, 2026
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
20 changes: 14 additions & 6 deletions pyfiles/Utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
import os
import subprocess

from osgeo import gdal, osr
from osgeo import gdal
from osgeo.gdalconst import GA_ReadOnly
import numpy as np

Expand Down Expand Up @@ -50,17 +50,25 @@ def initialize_output_raster_file(base_raster_file, output_raster_file, initial_
grid_initial_data = np.zeros((rows, cols), dtype=np.float32)
grid_initial_data[:] = float(initial_data)
else:
grid_initial_data = np.zeros((rows, cols), dtype=np.int)
grid_initial_data = np.zeros((rows, cols), dtype=np.int32)
grid_initial_data[:] = int(initial_data)

outband = outRaster.GetRasterBand(1)
outband.SetNoDataValue(NO_DATA_VALUE)
outband.WriteArray(grid_initial_data)

# set the projection of the tif file same as that of the base_raster file
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromWkt(base_raster.GetProjectionRef())
outRaster.SetProjection(outRasterSRS.ExportToWkt())
# Set projection from the source raster without reparsing WKT because
# some ArcGIS binary DEM projections can fail ImportFromWkt with corrupt data errors.
source_srs = base_raster.GetSpatialRef() if hasattr(base_raster, 'GetSpatialRef') else None
if source_srs:

@dtarb dtarb May 19, 2026

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not do the else option all the time?

@pkdash pkdash May 19, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added comments to the this code block to clarify the branching logic.

# GetSpatialRef is the preferred modern GDAL path and returns a spatial reference object that we export to WKT
outRaster.SetProjection(source_srs.ExportToWkt())
else:
# GetProjectionRef is retained only as a fallback for datasets/drivers where no spatial reference object is exposed
# Rely on dataset's raw projection string
source_wkt = base_raster.GetProjectionRef()
if source_wkt and source_wkt.strip():
outRaster.SetProjection(source_wkt)

outRaster = None

Expand Down
Loading