Skip to content
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.

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

UNA: allow the LULC raster to not have a nodata value defined #1631

Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
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
4 changes: 4 additions & 0 deletions HISTORY.rst
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,10 @@

Unreleased Changes
------------------
* Urban Nature Access
* The model now works as expected when the user provides an LULC raster
that does not have a nodata value defined.
https://github.com/natcap/invest/issues/1293
* Workbench
* Several small updates to the model input form UI to improve usability
and visual consistency (https://github.com/natcap/invest/issues/912)
Expand Down
45 changes: 39 additions & 6 deletions src/natcap/invest/urban_nature_access.py
Original file line number Diff line number Diff line change
Expand Up @@ -755,13 +755,12 @@ def execute(args):
squared_lulc_pixel_size = _square_off_pixels(args['lulc_raster_path'])

lulc_alignment_task = graph.add_task(
pygeoprocessing.warp_raster,
_warp_lulc,
kwargs={
'base_raster_path': args['lulc_raster_path'],
'target_pixel_size': squared_lulc_pixel_size,
'target_bb': target_bounding_box,
'target_raster_path': file_registry['aligned_lulc'],
'resample_method': 'near',
"source_lulc_path": args['lulc_raster_path'],
"target_lulc_path": file_registry['aligned_lulc'],
"target_pixel_size": squared_lulc_pixel_size,
"target_bounding_box": target_bounding_box,
},
target_path_list=[file_registry['aligned_lulc']],
task_name='Resample LULC to have square pixels'
Expand Down Expand Up @@ -2528,6 +2527,40 @@ def _create_mask(*raster_arrays):
_create_mask, target_mask_path, gdal.GDT_Byte, nodata_target=255)


def _warp_lulc(source_lulc_path, target_lulc_path, target_pixel_size,
target_bounding_box):
"""Warp a LULC raster and set a nodata if needed.

Args:
source_lulc_path (str): The path to a source LULC raster.
target_lulc_path (str): The path to the new LULC raster.
target_pixel_size (tuple): A 2-tuple of the target pixel size.
target_bounding_box (tuple): A 4-tuple of the target bounding box.

Returns:
``None``.
"""
source_raster_info = pygeoprocessing.get_raster_info(source_lulc_path)
target_nodata = source_raster_info['nodata'][0]
if target_nodata is None:
# Guarantee that our nodata cannot be represented by the datatype -
# select a nodata value that's out of range.
target_nodata = pygeoprocessing.choose_nodata(
source_raster_info['numpy_type']) + 1
Copy link
Member

Choose a reason for hiding this comment

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

I'm curious about the need to add 1 here. I mean, I understand from the comment why you've done this, but I wonder if you think choose_nodata should be updated to return the max + 1 instead of just the max? Would that make sense for other use cases, or could it cause problems?

Copy link
Member Author

Choose a reason for hiding this comment

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

Great question! Since we can assume that every pixel in the dataset is valid (because target_nodata is None), the safest possible nodata value is one that isn't even in the range of values possible with the given datatype. For a float32 raster, it's pretty unlikely that we would have a dataset that would have every discrete value that is allowed by the float32 spec, but it could still happen. Since the nodata value is stored as a piece of metadata separate from the pixel values themselves, we can avoid this possibility completely by using a nodata value that is completely outside of the range of possibilities of the datatype of the pixel values themselves.

If we were to change the datatype of choose_nodata, you are right that that would unfortunately break things. In raster_map, the nodata values of the incoming raster are transformed to the nodata values of the target raster, where the target nodata has been selected by choose_nodata. If choose_nodata were to return a value of the (dtype's max value)+1, then we'd get an overflow in the pixel values, or we bump up the output raster's datatype to accommodate the extra needed precision, neither of which we want!

So this solution as implemented gets around the situation by only handling this one case, where the lulc doesn't have a nodata value defined, and we just set one (to avoid breaking other things later on in the model) that can never happen in this datatype.


pygeoprocessing.warp_raster(
source_lulc_path, target_pixel_size, target_lulc_path,
'near', target_bb=target_bounding_box,
target_projection_wkt=source_raster_info['projection_wkt'])

# if there is no defined nodata, set a default value
raster = gdal.OpenEx(target_lulc_path, gdal.GA_Update)
band = raster.GetRasterBand(1)
band.SetNoDataValue(target_nodata)
Copy link
Member

Choose a reason for hiding this comment

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

Forgive me if I'm just misreading something, but could this SetNoDataValue unintentionally overwrite an existing nodata value that was defined already?

Copy link
Member Author

Choose a reason for hiding this comment

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

You know, it shouldn't be an issue because target_nodata would just be set to the same value it was at before, but this isn't really a clean approach. I just updated the function in 94f7f00 to only set the nodata value if one isn't already defined.

band = None
raster = None


def _mask_raster(source_raster_path, mask_raster_path, target_raster_path):
"""Convert pixels to nodata given an existing mask raster.

Expand Down
15 changes: 15 additions & 0 deletions tests/test_urban_nature_access.py
Original file line number Diff line number Diff line change
Expand Up @@ -358,6 +358,21 @@ def test_core_model(self):
self.assertAlmostEqual(numpy.min(valid_pixels), 1171.7352294921875)
self.assertAlmostEqual(numpy.max(valid_pixels), 11898.0712890625)

def test_no_lulc_nodata(self):
"""UNA: verify behavior when the LULC has no nodata value."""
from natcap.invest import urban_nature_access

args = _build_model_args(self.workspace_dir)
args['search_radius_mode'] = urban_nature_access.RADIUS_OPT_UNIFORM
args['search_radius'] = 100

raster = gdal.OpenEx(args['lulc_raster_path'], gdal.OF_RASTER)
band = raster.GetRasterBand(1)
band.DeleteNoDataValue()
band = None
raster = None
urban_nature_access.execute(args)

def test_split_urban_nature(self):
from natcap.invest import urban_nature_access

Expand Down
Loading