Align output NDR rasters to input DEM

[gtmaskall]

The output rasters from NDR should align with the input DEM coordinates. They currently do not. This, it seems, is because the model aligns to an intersection of the input raster and watershed bounding boxes.

The online documentation already states:

Raster inputs may have different cell sizes, and they will be resampled to match the cell size of the DEM. Therefore, all model results will have the same cell size as the DEM.

I propose, further, that the most consistent, desirable intended behaviour is that rasters are resampled onto the subset of DEM coordinates within whatever the final bounding box is. We're already treating (some) DEM spatial properties as somehow special.

Another strong reason to motivate this change is a use case wherein a user is processing two (or more) non-conterminous regions. It is entirely plausible, and more computationally efficient, to process the two regions separately. As the current resampling stands, this will all but guarantee the two sets of results lie on non-commensurable coordinates. That is, the multiple sets of coordinates will differ by an offset that is not a multiple of the pixel spacing.

Another justification is the case when the user has already done the work of aligning their input rasters. The current approach all but guarantees their outputs will lie on some offset coordinates. This seems ... "rude" ;-)

Further motivation centres on the subsequent use and analysis of the model outputs. The user may have other data, similarly aligned, they want to merge with the NDR output. They've now got additional work to do to cast the output back onto their preferred coordinates. It is true that the model provides aligned_ versions of the input rasters, but those inputs (required for driving the NDR model) are very likely not the only data the user has.

The user must ensure all input rasters are already in the same projected CRS. So, as with any software, there are already some expectations for data preparation on the user.

In short, the current approach, whilst attempting to resolve a "many to one" mapping of coordinates, actually inadvertently creates a many to many mapping when combining outputs across a number of runs (watersheds), and does this in a hard to predict and reconcile way, and also breaks what should be a "one to one" mapping of coordinates (in the case where all input rasters are aligned). Respecting the input DEM coordinates for the outputs will preserve the special case of "one to one" mapping for already aligned rasters, preserve the mapping across multiple runs (assuming the user has supplied consistent DEMs), and preserve the mapping a user might reasonably expect between their other raster data (that they will very likely have already aligned to the DEM) and the NDR model outputs.

[davemfish]

Well said @gtmaskall , I agree with these points.

We should take a look at all calls to pygeoprocessing.align_and_resize... across invest, not just NDR.