DEM Projections

[ecpease]

I've been delineating watersheds at my job for two years now using ArcGIS but would like to switch to Pysheds. Our method involves projecting the raw DEM to an Albers Equal Area projection from geographical coordinate system first, then continuing with fill, flow direction, etc. This is "the way we do it here" but I know Pysheds doesn't involve projections before raster processing. From what I've done and read, it doesn't seem to matter much whether projection occurs before or after. Would you mind describing your reasoning for not projecting so that I can understand the process better? Does it actually matter or make a difference? I just really want to understand the process as thoroughly as I can. Thank you in advance!

[mdbartos]

Hi Emily,

The main difference that I can see is that projecting will affect the distances between adjacent cells, and therefore will affect the slope computation. While the difference is small, I think projecting first is technically correct (unless the data provider has already accounted for this, which I think HydroSHEDs does).

However, projecting beforehand comes with a performance penalty. If the original raster is defined on a regular grid, projecting from geographic coordinates will result in an irregularly-spaced set of points. This means that a lot of the array operations I used to speed things up will no longer work.

There are two ways I can think of to address this problem:

1. Project, and then keep track of all the distances between adjacent cells when computing the slopes.

2. Project, and then resample the projected points to a regularly-spaced grid in the new coordinate system.

---

I should hopefully have a code sample using solution (1) posted soon.

[mdbartos]

Hi Emily,

I've added the ability compute flow directions under a given projection. You can apply the projection by using the as_crs keyword argument in the call to grid.flowdir. For example:

new_crs = pyproj.Proj('+init=epsg:3083')
grid.flowdir('inflated_dem', out_name='dir', dirmap=dirmap, as_crs=new_crs)

The following notebook compares the projected and unprojected cases:

https://github.com/mdbartos/pysheds/blob/master/examples/flowdir_vs_projection.ipynb

While there are some differences, they appear to be pretty small.