NDR raster output alignment

Hi,

I’m running the NDR model and my DEM, Land Use and Nutrient Runoff Proxy rasters all have the same projection, with their cells lining up perfectly with a regular grid, but it seems that the main retention output aligns the extent to that of the Watersheds vector (in the same projection as the rasters). As such, the output raster is not aligned with the regular OS grid, and is unsuitable for further analysis in that state. Is there a way to change the output extent so that it aligns with a regular grid, i.e. 100 x 100m raster cells sit neatly within a 1km grid square rather than are off by several metres?

I’ve tried using one inland watershed instead of the whole vector, which results in the cells lining up - well, they’re 0.001m to the left and down but good enough.

I’ve also tried using the single most southwestern watershed, which results in the exact same irregular alignment as when using the whole watersheds vector (below). This makes me think it’s possibly to do with how the vector overlaps the raster layers in this southwestern watershed, not necessarily where the watersheds vector origin lies.

I’ve attached the logfile of the southwestern watershed (the one that shares the same origin as the entire watersheds vector). I’d appreciate any help offered, as I’m not really sure what to try next beyond wading into the python script!

Thanks :slight_smile:

InVEST-Nutrient-Delivery-Ratio-log-2023-10-02–16_09_44.txt (22.5 KB)

1 Like

Hi @emilyu , thanks for your illustration of the problem, this is an interesting case.

The NDR model is trying to be helpful by handling the case where none of the input rasters are aligned with each other, but in your case they are already well-aligned, so it’s not being very helpful. Specifically, the model does this:

pygeoprocessing.align_and_resize_raster_stack(
    base_raster_list,
    aligned_raster_list,
    ['near']*len(base_raster_list),
    dem_info['pixel_size'],
    'intersection',
    'base_vector_path_list'=[args['watersheds_path']])

This operation is taking all the input raster datasets, and the watersheds vector, and taking the intersection of all those layers’ bounding boxes. Then the output rasters are aligned to that “bounding-box-intersection” (specifically the upper-left corner of it, I think). So I think that’s why you see slightly different behavior when you tried those different watersheds.

It’s interesting to hear that incorporating the bounding box of the watersheds layer was actually unhelpful in your case. We could consider modifying the model to force the output grid to always align to one of the specific input rasters (while still being completely contained within that “bounding-box-intersection”). Unfortunately there’s no way for a user to choose that option right now, without modifying the source code.

However, it’s worth noting that the output workspace includes “aligned” copies of all your input datasets. So if you are wanting to do some post-processing where you use your input rasters together with the output rasters, you could consider using those aligned versions of the inputs. You can find them in intermediate_outputs/cache_dir/ with filenames like aligned_*.tif. Does that meet your needs?

1 Like

@dave, for what it’s worth, I have also had occasions where the change in alignment of InVEST outputs is a problem, because I need the model output pixels to align exactly with non-InVEST rasters. This doesn’t happen often to me, but it’s unfortunately hard to remedy in post-processing, and any re-aligning in post-processing might lead to errors.

So I would also support being able to have one layer (like the LULC, or DEM, perhaps whichever one is currently used for reprojection, or user-selectable) whose alignment is unchanged, and used for snapping of the others. But then I wonder why, in the cases where we have already pre-processed our inputs to be in the same projection, and aligned, InVEST does it again no matter what. At the moment, there’s probably no way for the model to know that they are aligned, so I’m guessing that it does the most conservative thing?

~ Stacie

1 Like

I agree this makes a lot of sense for the cases where we are also defining the target pixel size from one of these inputs, which we probably are in most/all cases.

I guess it’s just easier to do the alignment than to check whether they are already aligned. And I’m guessing that somewhere deeper down GDAL is optimized such that it won’t do any expensive warping if the projection, pixel size, and bounding box are all the same going in as they will be coming out.

But in some cases, like NDR, align_and_resize is also using the bounding box of the watersheds vector as a way to clip the input rasters to avoid some extra computation. So in theory, that’s useful even when the rasters are already aligned. Anyhow, we can and should see if there’s a better way to do this that includes snapping to one of the input rasters.

1 Like

Thanks for your comments, @dave and @swolny :slight_smile:

I found a way round the problem, for now at least. I’ve opted for a quick and dirty method of adding a 1km square “watershed” to my watersheds vector at a more useful southwest origin - it doesn’t overlap the input rasters so it’s not calculating an export value there, but it means everything is lining up as it should! It works for my data but might not necessarily for everyone’s.

In the future, it would be great to be able to choose how the export raster aligns. For me, at least, this output forms part of a bigger analysis, and so it’s very important that the export aligns exactly with other datasets.

Thanks again for your swift comments!

2 Likes

Ah ha!! Y’all have been talking about this issue! :slight_smile:

This has bitten me, also. Having read this thread and checked for any GH issues, I’ve taken the liberty of raising one. I hope that helps support the case for change.

As for myself, my workflow broadly comprises:

  1. Start with vector of desired watershed target area
  2. Pull DEM to cover that area
  3. Prepare other required inputs and align to the DEM raster
  4. Run ndr.execute
  5. Analyze output in the context of other (often) raster data - so wanting coordinates to align is quite common

I get the feeling, from this forum, that my coming from the direction of python, the STAC API, xarray etc. is not super common. For those users of xarray I recommend its interp_like method to interpolate the NDR output(s) onto your expected DEM coordinates.

2 Likes

This topic was automatically closed 90 days after the last reply. New replies are no longer allowed.