Accumulation is wrongly calculated as all 1s.

[bosmanoglu]

First, I would like to thank you very much for developing this tool. It's well written and documented allowing me to integrate in one of our projects.

I am working on calculating HAND for a given area with detected flood extent. It's a fairly large area so I am separating each contiguous flood extent as one object and processing a rectangular area encompassing each object individually. Following your examples, I am using the following code to process each of these:

for l  in labels_of_labeled_flood_mask:
  #some code to clip dem to rectangular area given the labeled flood mask
  #this gives me variables: dem, dem_proj4, dem_gT
  #dem = https://www.dropbox.com/s/zkib6jlfejuc1vi/dem.npy?dl=0
  #dem_proj4 = Proj('+proj=longlat +datum=WGS84 +no_defs', preserve_units=True)
  #dem_gT = (90.9230473, 0.0002873612583231426, 0.0, 25.4224585, 0.0, -0.00028736127733026517)
  # clip_gT()=(90.99603705961407, 0.0002873612583231426, 0.0, 25.369584024971232, 0.0, -0.00028736127733026517)
  #cip_gT() calculates clipped geoTransform (gT) based on image coordinates, output above.
  dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
  affine=Affine.from_gdal(*clip_gT(dem_gT, min0, max0, min1, max1, method='image'))
  mask=np.ones(dem.shape, dtype=np.bool)
  grid=Grid(shape=dem.shape,affine=affine, crs=dem_proj4, mask=mask)
  grid.add_gridded_data(dem, data_name='dem',affine=affine, crs=dem_proj4, mask=mask)
  grid.fill_depressions('dem', out_name='flooded_dem')
  grid.resolve_flats('flooded_dem', out_name='inflated_dem')
  grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap)
  grid.accumulation(data='inflated_dem', dirmap=dirmap, out_name='acc')
  grid.compute_hand('dir', 'inflated_dem', grid.acc >0, out_name='hand')
  hand_output=grid.view('hand')

Here is a look at the DEM:

What I notice is that the accumulation result is all 1, as seen below:

ipdb> (grid.view('acc')==1).all()
Raster(True)

The flow_dir looks reasonable though:

Because accumulation was all 1's I set the grid.acc >0 in compute_hand calculation, despite noticing that in the test_grid.py you have it set as grid.acc > 100. As a result HAND is all 0's except a single row/column of nan's surrounding it.

ipdb> grid.view('hand')
Raster([[nan, nan, nan, ..., nan, nan, nan],
        [nan,  0.,  0., ...,  0.,  0., nan],
        [nan,  0.,  0., ...,  0.,  0., nan],
        ...,
        [nan,  0.,  0., ...,  0.,  0., nan],
        [nan,  0.,  0., ...,  0.,  0., nan],
        [nan, nan, nan, ..., nan, nan, nan]])

It's likely that I did something wrong, and there is nothing wrong with pysheds. Could you please take a look and let me know which is true? If my method/input is correct but there is a glitch in pysheds, if you can give me pointers to look into, I'm happy to check and let you know what I find.

[mdbartos]

The accumulation function takes the flow direction file as input. Therefore:

grid.accumulation(data='inflated_dem', dirmap=dirmap, out_name='acc')

Should be:

grid.accumulation(data='dir', dirmap=dirmap, out_name='acc')

[bosmanoglu]

Yes, it looks like I made a copy paste error. With your suggested correction we were able to get results.

This ticket can be closed.