Skip to content

Commit

Permalink
Fix Issue #121: force layer to disk before saving.
Browse files Browse the repository at this point in the history
  • Loading branch information
@mdales
mdales committed Aug 20, 2024
1 parent 54260c1 commit d159f27
Show file tree
Hide file tree
Showing 2 changed files with 113 additions and 129 deletions.
240 changes: 112 additions & 128 deletions methods/inputs/generate_slope.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,138 +111,122 @@ def generate_slope(input_elevation_directory: str, output_slope_directory: str):
continue

with tempfile.TemporaryDirectory() as tmpdir:
elevation = RasterLayer.layer_from_file(elev_path)

logging.info("Area of elevation tile %a", elevation.area)
_easting, _northing, lower_code, lower_letter = utm.from_latlon(
elevation.area.bottom, elevation.area.left
)
_easting, _northing, upper_code, upper_letter = utm.from_latlon(
elevation.area.top, elevation.area.right
)

# FAST PATH -- with only one UTM zone the reprojection back has no issues
if lower_code == upper_code and lower_letter == upper_letter:
actual_utm_code = lower_code
warp(
actual_utm_code,
elev_path,
elevation.pixel_scale.xstep,
elevation.pixel_scale.ystep,
out_path,
with RasterLayer.layer_from_file(elev_path) as elevation:

logging.info("Area of elevation tile %a", elevation.area)
_easting, _northing, lower_code, lower_letter = utm.from_latlon(
elevation.area.bottom, elevation.area.left
)
else:
# SLOW PATH -- in the slow path, we have to break the elevation raster into
# UTM sections and do the above to each before reprojecting back and recombining

# To capture the results here for later inspection just override the tmpdir variable
for actual_utm_code in range(lower_code, upper_code + 1):
for utm_letter in crange(lower_letter, upper_letter):
logging.debug("UTM(%s,%s)", actual_utm_code, utm_letter)

# Note: we go a little bit around the UTM tiles and will crop them down to size later
# this is to remove some aliasing effects.
bbox = bounding_box_of_utm(actual_utm_code, utm_letter, UTM_EXPANSION_DEGREES)

# Crop the elevation tile to a UTM zone
utm_layer = RasterLayer.empty_raster_layer_like(
elevation, area=bbox
)
utm_id = f"{actual_utm_code}-{utm_letter}-{elevation_path}"
utm_clip_path = os.path.join(tmpdir, utm_id)
intersection = RasterLayer.find_intersection(
[elevation, utm_layer]
)
result = RasterLayer.empty_raster_layer(
intersection,
elevation.pixel_scale,
elevation.datatype,
utm_clip_path,
elevation.projection,
)
result.set_window_for_intersection(intersection)
elevation.set_window_for_intersection(intersection)
elevation.save(result)

# Flush elevation utm clip to disk
del result

# Now warp into UTM, calculate slopes, and warp back
slope_out_path = os.path.join(tmpdir, "out-slope-" + utm_id)
warp(
actual_utm_code,
utm_clip_path,
elevation.pixel_scale.xstep,
elevation.pixel_scale.ystep,
slope_out_path,
)

# We now recrop the out-slope back to the bounding box we assumed at the start
bbox_no_expand = bounding_box_of_utm(
actual_utm_code, utm_letter, 0.0
)
slope_tif = RasterLayer.layer_from_file(slope_out_path)
grid = RasterLayer.empty_raster_layer_like(
slope_tif, area=bbox_no_expand
)
output_final = f"final-slope-{actual_utm_code}-{utm_letter}-{elevation_path}"
final_path = os.path.join(tmpdir, output_final)
logging.debug("Slope underlying %s", slope_tif._underlying_area) # pylint: disable=W0212
logging.debug("Grid underling %s", grid._underlying_area) # pylint: disable=W0212
try:
intersection = RasterLayer.find_intersection([slope_tif, grid])
except ValueError:
logging.debug(
"UTM (%s, %s) didn't intersect actual area %s",
actual_utm_code,
utm_letter,
grid._underlying_area # pylint: disable=W0212
)
continue
slope_tif.set_window_for_intersection(intersection)
final = RasterLayer.empty_raster_layer(
intersection,
slope_tif.pixel_scale,
slope_tif.datatype,
final_path,
slope_tif.projection,
)
logging.debug("Final underlying %s", final._underlying_area) # pylint: disable=W0212
final.set_window_for_intersection(intersection)
slope_tif.save(final)

# Flush
del final

# Now to recombine the UTM gridded slopes into the slope tile
slopes = glob("final-slope-*", root_dir=tmpdir)
assert len(slopes) > 0

# This sets the order a little better for the union of the layers
slopes.sort()
slopes.reverse()

logging.info("Render order %s", slopes)

combined = GroupLayer(
[
RasterLayer.layer_from_file(os.path.join(tmpdir, filename))
for filename in slopes
]
_easting, _northing, upper_code, upper_letter = utm.from_latlon(
elevation.area.top, elevation.area.right
)

elevation = RasterLayer.layer_from_file(elev_path)
intersection = RasterLayer.find_intersection([elevation, combined])
combined.set_window_for_intersection(intersection)
elevation.set_window_for_intersection(intersection)

assembled_path = os.path.join(tmpdir, "patched.tif")
result = RasterLayer.empty_raster_layer_like(
elevation, filename=assembled_path
)
combined.save(result)
# FAST PATH -- with only one UTM zone the reprojection back has no issues
if lower_code == upper_code and lower_letter == upper_letter:
actual_utm_code = lower_code
warp(
actual_utm_code,
elev_path,
elevation.pixel_scale.xstep,
elevation.pixel_scale.ystep,
out_path,
)
else:
# SLOW PATH -- in the slow path, we have to break the elevation raster into
# UTM sections and do the above to each before reprojecting back and recombining

# To capture the results here for later inspection just override the tmpdir variable
for actual_utm_code in range(lower_code, upper_code + 1):
for utm_letter in crange(lower_letter, upper_letter):
logging.debug("UTM(%s,%s)", actual_utm_code, utm_letter)

# Note: we go a little bit around the UTM tiles and will crop them down to size later
# this is to remove some aliasing effects.
bbox = bounding_box_of_utm(actual_utm_code, utm_letter, UTM_EXPANSION_DEGREES)

# Crop the elevation tile to a UTM zone
with RasterLayer.empty_raster_layer_like(elevation, area=bbox) as utm_layer:
utm_id = f"{actual_utm_code}-{utm_letter}-{elevation_path}"
utm_clip_path = os.path.join(tmpdir, utm_id)
intersection = RasterLayer.find_intersection(
[elevation, utm_layer]
)
with RasterLayer.empty_raster_layer(
intersection,
elevation.pixel_scale,
elevation.datatype,
utm_clip_path,
elevation.projection,
) as result:
result.set_window_for_intersection(intersection)
elevation.set_window_for_intersection(intersection)
elevation.save(result)

# Now warp into UTM, calculate slopes, and warp back
slope_out_path = os.path.join(tmpdir, "out-slope-" + utm_id)
warp(
actual_utm_code,
utm_clip_path,
elevation.pixel_scale.xstep,
elevation.pixel_scale.ystep,
slope_out_path,
)

shutil.move(assembled_path, out_path)
# We now recrop the out-slope back to the bounding box we assumed at the start
bbox_no_expand = bounding_box_of_utm(
actual_utm_code, utm_letter, 0.0
)
with RasterLayer.layer_from_file(slope_out_path) as slope_tif:
with RasterLayer.empty_raster_layer_like(slope_tif, area=bbox_no_expand) as grid:
output_final = f"final-slope-{actual_utm_code}-{utm_letter}-{elevation_path}"
final_path = os.path.join(tmpdir, output_final)
logging.debug("Slope underlying %s", slope_tif._underlying_area) # pylint: disable=W0212
logging.debug("Grid underling %s", grid._underlying_area) # pylint: disable=W0212
try:
intersection = RasterLayer.find_intersection([slope_tif, grid])
except ValueError:
logging.debug(
"UTM (%s, %s) didn't intersect actual area %s",
actual_utm_code,
utm_letter,
grid._underlying_area # pylint: disable=W0212
)
continue
slope_tif.set_window_for_intersection(intersection)
with RasterLayer.empty_raster_layer(
intersection,
slope_tif.pixel_scale,
slope_tif.datatype,
final_path,
slope_tif.projection,
) as final:
logging.debug("Final underlying %s", final._underlying_area) # pylint: disable=W0212
final.set_window_for_intersection(intersection)
slope_tif.save(final)

# Now to recombine the UTM gridded slopes into the slope tile
slopes = glob("final-slope-*", root_dir=tmpdir)
assert len(slopes) > 0

# This sets the order a little better for the union of the layers
slopes.sort()
slopes.reverse()

logging.info("Render order %s", slopes)

with GroupLayer.layer_from_files([os.path.join(tmpdir, filename) for filename in slopes]) as combined:
with RasterLayer.layer_from_file(elev_path) as elevation:
intersection = RasterLayer.find_intersection([elevation, combined])
combined.set_window_for_intersection(intersection)
elevation.set_window_for_intersection(intersection)

assembled_path = os.path.join(tmpdir, "patched.tif")
with RasterLayer.empty_raster_layer_like(
elevation, filename=assembled_path
) as result:
combined.save(result)

shutil.move(assembled_path, out_path)


def main() -> None:
Expand Down
2 changes: 1 addition & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ scipy
numba
matplotlib
geojson
git+https://github.com/quantifyearth/yirgacheffe@bd2e91c773a414f66340ebb8c13044a1b1a6045f
git+https://github.com/quantifyearth/yirgacheffe@cc89b4d8a0e97c1a11b730cd688a58b680023336
git+https://github.com/carboncredits/biomass-recovery@9e54f80832a7eca915ebd13b03df9db2a08aee9d

# developement
Expand Down

0 comments on commit d159f27

Please sign in to comment.