r-barnes / richdem

High-performance Terrain and Hydrology Analysis
GNU General Public License v3.0
266 stars 68 forks source link

Error in FlowAccumulation #7

Closed theobarnhart-USGS closed 6 years ago

theobarnhart-USGS commented 6 years ago

Hi Rich,

I get the following error in richdem.FlowAccumulation with both the 'D8' and 'Dinf' methods. Any thoughts on how to resolve it are appreciated.

Thank you!

In [7]: dem = rd.LoadGDAL('./cedar2m/hdr.adf')

In [8]: rd.FillDepressions(dem,in_place=True)

A Priority-Flood (Zhou2016 version)
C Zhou, G., Sun, Z., Fu, S., 2016. An efficient variant of the Priority-Flood algorithm for filling depressions in raster digital el
evation models. Computers & Geosciences 90, Part A, 87 – 96. doi:http://dx.doi.org/10.1016/j.cageo.2016.02.021

t Zhou2016 wall-time = 169.659 s

In [9]: accumD8 = rd.FlowAccumulation(dem,method='D8')

A Generic Flow Accumulation Algorithm

A O'Callaghan (1984)/Marks (1984) Flow Accumulation (aka D8)
C O'Callaghan, J.F., Mark, D.M., 1984. The Extraction of Drainage Networks from Digital Elevation Data. Computer vision, graphics, a
nd image processing 28, 323--344.

[                                                  ] (0% - infs - 1 threads)--------------------------------------------------------
-------------------
IndexError                                Traceback (most recent call last)
<ipython-input-9-402d8b713134> in <module>()
----> 1 accumD8 = rd.FlowAccumulation(dem,method='D8')

~/.conda/envs/py36/lib/python3.6/site-packages/richdem/__init__.py in FlowAccumulation(dem, method, exponent, weights, in_place)
    422
    423   if method in facc_methods:
--> 424     facc_methods[method](dem.wrap(),accumw)
    425   elif method in facc_methods_exponent:
    426     if exponent is None:

IndexError: vector::_M_range_check
theobarnhart-USGS commented 6 years ago

The same error also occurs with the Rho8 method as well.

theobarnhart-USGS commented 6 years ago

Hi Richard,

I just wanted to check in on this issue, I updated richDEM today and the above error still persists. Any idea about what may be causing it? Is it something weird in my input data?

Thanks!

r-barnes commented 6 years ago

Sorry, I've just seen this. I haven't been able to reproduce the problem on my machine. Could you verify that it is persisting?

theobarnhart-USGS commented 6 years ago

Hi, I just re-ran this and it looks like the issue is still persisting.

(py36) [tbarnhart@n3-134 cedar2m]$ pip install --upgrade richdem
Requirement already up-to-date: richdem in /home/tbarnhart/.conda/envs/py36/lib/python3.6/site-packages
Collecting numpy<2,>=1.7; python_version > "3.4" or python_version < "3.0" (from richdem)
  Downloading https://files.pythonhosted.org/packages/71/90/ca61e203e0080a8cef7ac21eca199829fa8d997f7c4da3e985b49d0a107d/numpy-1.14.3-cp36-cp36m-manylinux1_x86_64.whl (12.2MB)
    100% |████████████████████████████████| 12.2MB 41kB/s
Installing collected packages: numpy
  Found existing installation: numpy 1.14.2
    Uninstalling numpy-1.14.2:
      Successfully uninstalled numpy-1.14.2
Successfully installed numpy-1.14.3
You are using pip version 9.0.3, however version 10.0.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.
(py36) [tbarnhart@n3-134 cedar2m]$ ipython
Python 3.6.2 |Continuum Analytics, Inc.| (default, Jul 20 2017, 13:51:32)
Type 'copyright', 'credits' or 'license' for more information
IPython 6.2.1 -- An enhanced Interactive Python. Type '?' for help.

Caching the list of root modules, please wait!
(This will only be done once - type '%rehashx' to reset cache!)

In [1]: import richdem as rd

In [2]: dem = rd.LoadGDAL('./hdr.adf'
   ...: )

In [3]: rd.FillDepressions(dem, in_place=True)

A Priority-Flood (Zhou2016 version)
C Zhou, G., Sun, Z., Fu, S., 2016. An efficient variant of the Priority-Flood algorithm for filling depressions in raster digital elevation models. Computers & Geosciences 90, Part A, 87 – 96. doi:http://dx.doi.org/10.1016/j.cageo.2016.02.021

t Zhou2016 wall-time = 254.165 s

In [4]: accumD8 = rd.FlowAccumulation(dem,method='D8')

A Generic Flow Accumulation Algorithm

A O'Callaghan (1984)/Marks (1984) Flow Accumulation (aka D8)
C O'Callaghan, J.F., Mark, D.M., 1984. The Extraction of Drainage Networks from Digital Elevation Data. Computer vision, graphics, and image processing 28, 323--344.

[                                                  ] (0% - infs - 1 threads)---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-4-402d8b713134> in <module>()
----> 1 accumD8 = rd.FlowAccumulation(dem,method='D8')

~/.conda/envs/py36/lib/python3.6/site-packages/richdem/__init__.py in FlowAccumulation(dem, method, exponent, weights, in_place)
    453
    454   if method in facc_methods:
--> 455     facc_methods[method](dem.wrap(),accumw)
    456   elif method in facc_methods_exponent:
    457     if exponent is None:

IndexError: vector::_M_range_check
r-barnes commented 6 years ago

Could you send me the data file you're working with? Although the algorithms should work with any data type, there are a few I haven't done extensive testing on. I have a feeling that might be the case here.

theobarnhart-USGS commented 6 years ago

Here is a link to the data file: ftp://dnrftp.dnr.ne.gov/Pub/data/dems/LiDAR_2011/Cedar2m.zip

r-barnes commented 6 years ago

Thanks. I can reproduce the issue using the following code:

import richdem as rd
dem = rd.LoadGDAL('./cedar2m/hdr.adf')
rd.FillDepressions(dem, in_place=True)
accumD8 = rd.FlowAccumulation(dem,method='D8')

I'm looking into it.

r-barnes commented 6 years ago

The following C++ code reproduces the issue, making it easier to debug than in Python:

#include "richdem/common/Array2D.hpp"
using namespace richdem;

int main(int argc, char **argv){
  Array2D<float> dem(argv[1]);
  Array2D<double> accum(dem,1);
  FA_D8(dem,accum);
}

It can be compiled with:

GDAL_LIBS="`gdal-config --libs`"
GDAL_CFLAGS="`gdal-config --cflags` -DUSEGDAL"
#RICHDEM_GIT_HASH=`git rev-parse HEAD`
#RICHDEM_COMPILE_TIME=`date -u +'%Y-%m-%d %H:%M:%S UTC'`
CXXFLAGS="$GDAL_CFLAGS --std=c++11 -g -O0 -Wall -Wno-unknown-pragmas -Irichdem/include"
LIBS="$GDAL_LIBS"

g++ $CXXFLAGS -o temp.exe temp.cpp $LIBS -Wno-sign-compare
r-barnes commented 6 years ago

Confirmed that the bug does not appear on the 2418x1636 Float32 beauford dataset, but does appear on the 19742x28685 Float32 cedar dataset.

r-barnes commented 6 years ago

The bug turned out to arise from an integer overflow error in resizing the flow proportions array. That is, the input data was larger than the flow accumulation code expected, though still well within the operable range of the software.

This bug is fixed as of 9b523dadd43816dbbc08a20873bed057b1f97048 on the dev branch.

@theobarnhart-USGS: Are you up for testing from dev? (If you're on Linux/Mac, this requires cloning the repo and running python3 setup.py install --user in the wrappers/pyrichdem directory.) Alternatively, I can wrap up some stuff I'm working on and release an updated Python version which you can get via pip.

Also: I notice that you're filling depressions without an epsilon gradient or flat resolution. The result may not be what you want.

theobarnhart-USGS commented 6 years ago

I think I should be able to test that from dev. I can probably give it a shot tomorrow. Thank you for investigating that bug!

On Mon, May 7, 2018 at 5:37 PM, Richard Barnes notifications@github.com wrote:

The bug turned out to arise from an integer overflow error in resizing the flow proportions array. That is, the input data was larger than the flow accumulation code expected, though still well within the operable range of the software.

This bug is fixed as of 9b523da https://github.com/r-barnes/richdem/commit/9b523dadd43816dbbc08a20873bed057b1f97048 on the dev branch.

@theobarnhart-USGS https://github.com/theobarnhart-USGS: Are you up for testing from dev? (If you're on Linux/Mac, this requires cloning the repo and running python3 setup.py install --user in the wrappers/pyrichdem directory.) Alternatively, I can wrap up some stuff I'm working on and release an updated Python version which you can get via pip.

Also: I notice that you're filling depressions without an epsilon gradient https://richdem.readthedocs.io/en/latest/depression_filling.html#epsilon-filling or flat resolution https://richdem.readthedocs.io/en/latest/flat_resolution.html. The result may not be what you want.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/r-barnes/richdem/issues/7#issuecomment-387238816, or mute the thread https://github.com/notifications/unsubscribe-auth/Ac3kGksWGkssdOvhg35fWdxHywYkjwvPks5twNq4gaJpZM4RtQtj .

-- Theodore Barnhart tbarnhart@usgs.gov Physical Scientist WY-MT Water Science Center U.S. Geological Survey Helena, MT

theobarnhart-USGS commented 6 years ago

I can't really tell if I've built the dev version correctly, but here is the error that it throws now:

In [4]: accumD8 = rd.FlowAccumulation(dem,method='D8')

A Generic Flow Accumulation Algorithm

A O'Callaghan (1984)/Marks (1984) Flow Accumulation (aka D8)
C O'Callaghan, J.F., Mark, D.M., 1984. The Extraction of Drainage Networks from Digital Elevation Data. Computer vision, graphics, and image processing 28, 323--344.

[                                                  ] (0% - infs - 1 threads)---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-4-402d8b713134> in <module>()
----> 1 accumD8 = rd.FlowAccumulation(dem,method='D8')

~/.local/lib/python3.6/site-packages/richdem-0.2.0-py3.6-linux-x86_64.egg/richdem/__init__.py in FlowAccumulation(dem, method, exponent, weights, in_place)
    453
    454   if method in facc_methods:
--> 455     facc_methods[method](dem.wrap(),accumw)
    456   elif method in facc_methods_exponent:
    457     if exponent is None:

IndexError: vector::_M_range_check: __n (which is 801726138) >= this->size() (which is 801726134)
r-barnes commented 6 years ago

@theobarnhart-USGS: I'm sorry about the difficulties here.

I've released a new version of the code to PyPI. It should be installable with pip install richdem and, hopefully, work. Could you test again?

I've also overhauled the release system to make it easier to push out bug fixes to issues like this.

theobarnhart-USGS commented 6 years ago

I'll give it a shot, thank you Rich!

On Mon, May 28, 2018 at 4:20 PM, Richard Barnes notifications@github.com wrote:

@theobarnhart-USGS https://github.com/theobarnhart-USGS: I'm sorry about the difficulties here.

I've released a new version of the code to PyPI. It should be installable with pip install richdem and, hopefully, work. Could you test again?

I've also overhauled the release system to make it easier to push out bug fixes to issues like this.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/r-barnes/richdem/issues/7#issuecomment-392611224, or mute the thread https://github.com/notifications/unsubscribe-auth/Ac3kGuKwIJVehKcJzxMEWOUPA8RZI7m3ks5t3HgrgaJpZM4RtQtj .

-- Theodore Barnhart, Ph.D. tbarnhart@usgs.gov Physical Scientist WY-MT Water Science Center U.S. Geological Survey Helena, MT

theobarnhart-USGS commented 6 years ago

version 0.3.3 appears to have resolved this. Thank you @r-barnes!

In [17]: dem = rd.LoadGDAL('./data/NE/cedar2m/hdr.adf')

In [18]: accum = rd.FlowAccumulation(dem,method='D8')

A O'Callaghan (1984)/Marks (1984) D8/D4 Flow Accumulation
C O'Callaghan, J.F., Mark, D.M., 1984. The Extraction of Drainage Networks from Digital Elevation Data. Computer vision, graphics, and image processing 28, 323--344.

c topology = D8

A Generic Flow Accumulation Algorithm
p Creating dependencies array...

d Source cells found = 144528093
p Calculating flow accumulation...
t Wall-time       = 53.3961 s