Hi List,
I'm a grass newbie. I'm trying to route some water through a DEM and think r.terraflow is the right tool. It is crashing. Can someone point out what I'm doing wrong? I have not found much when searching for this problem on google or list archives. I think there are several issues occurring here. I'm having the same problems on grass64 (installs easily on OS X) and grass71 (on a linux machine).
I'm working in the following region:
g.region -p
projection: 99 (Stereographic)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: -657600
south: -3349500
west: -638000
east: 864700
nsres: 150
ewres: 150
rows: 17946
cols: 10018
cells: 179783028
1) r.mapcalc produces file of 0 size
When I run
r.mapcalc 'phi = (1000*9.8*bed + 917*9.8*(surf-bed))'
The file has 0 size. I now do the same calculation using nctools and then load the netcdf file, but I'd rather do this in grass if I could. FYI, bed is 589 MB and surf is 190 MB. When I calculate phi with nctools (ncap2 -O -v -s "phi=1E3*9.8*bed+917*9.8*thickness" in.nc phi.nc) the resulting file is 1.4 GB. Clearly the 0 size grass-produced phi is due to overflow int/float/double issues? But a) there is no warning/error and b) wrapping the variables above in float()'s doesn't seem to help. The file still has 0 size if I do a similar but simpler computation:
r.mapcalc 'phi = (bed + 0.917(surf-bed))'
2) r.terraflow crashes with 'Bad address' error.
r.terraflow --o --v elevation=phi fill=f direction=d sw=s acc=a tci=tci
The results of that command:
cell phi header compatible with region header
Elevation stored as FLOAT (4B)
WARNING: raster phi is of type CELL_TYPE --you should use r.terraflow.short
Region size is 17946 x 10018
STREAM temporary files in /home/user/data/grass/GR/PERMANENT/.tmp/host/26900.0 (THESE INTERMEDIATE STREAMS WILL NOT BE DELETED IN CASE OF ABNORMAL TERMINATION OF THE PROGRAM. TO SAVE SPACE PLEASE DELETE THESE FILES MANUALLY!)
Memory manager registering memory in MM_WARN_ON_MEMORY_EXCEEDED mode.
Reading data from <phi> to stream
</home/user/data/grass/GR/PERMANENT/.tmp/host/26900.0/STREAM_baaaaX>
100%
total elements=179783028, nodata elements=18914292
largest temporary files:
FILL: 9.38G (10067849568) [179783028 elements, 56B each]
FLOW: 11.99G (12869498880) [160868736 elements, 80B each]
Will need at least 23.97G (25738997760) space available in
/home/user/data/grass/GR/PERMANENT/.tmp/host/26900.0
MM warning: limit=314572800B. allocating 156610272B. limit exceeded by 3706B.
EMPQUEUEADAPTIVE: starting in-memory pqueue
EMPQUEUEADAPTIVE: available memory: 297.937MB
EMPQUEUEADAPTIVE: desired memory: 297.937MB
sz_stream: 270400 buf_arity: 200 mm_overhead: 8666496 mm_avail: 312409214.
EMPQUEUEADAPTIVE: memory overhead set to 8.26501MB
EMPQUEUEADAPTIVE: pqsize set to 37967839
EMPQUEUEADAPTIVE: starting in-memory pqueue
EMPQUEUEADAPTIVE: available memory: 297.163MB
EMPQUEUEADAPTIVE: desired memory: 297.163MB
sz_stream: 270400 buf_arity: 200 mm_overhead: 8666496 mm_avail: 311597910.
EMPQUEUEADAPTIVE: memory overhead set to 8.26501MB
EMPQUEUEADAPTIVE: pqsize set to 37866426
EMPQUEUEADAPTIVE: starting in-memory pqueue
EMPQUEUEADAPTIVE: available memory: 296.389MB
EMPQUEUEADAPTIVE: desired memory: 296.389MB
sz_stream: 270400 buf_arity: 200 mm_overhead: 8666496 mm_avail: 310786606.
EMPQUEUEADAPTIVE: memory overhead set to 8.26501MB
EMPQUEUEADAPTIVE: pqsize set to 37765013
MM warning: limit=314572800B. allocating 155798940B. limit exceeded by 3706B.
rusage/gettimeofday: Bad address
Thanks for any help,
-k.
On 2016-01-03 at 15:46, Anna Petrášová <kratochanna@gmail.com> wrote:
I don't have any experience with r.terraflow, but try using
r.watershed.
Crashes with different errors. I'll post about that later. I was using terraflow because it claims to be for large data (not that I think this data is particularly large...)
r.mapcalc 'phi = (1000*9.8*bed + 917*9.8*(surf-bed))'
no idea, have you tried testing it on a smaller region?
downsampling doesn't help.
r.mapcalc 'phi = bed*2' does produce a data file, full resolution or downsampled.
Oddly, "2*bed" works, but "2.5*bed" does not work.
Does r.mapcalc only do integer math? Not according to the manual. "r.info bed" shows the range is -5028 to 3678.
-k.
On Sun, Jan 3, 2016 at 11:55 PM, Ken Mankoff <mankoff@gmail.com> wrote:
On 2016-01-03 at 15:46, Anna Petrášová <kratochanna@gmail.com> wrote:
I don't have any experience with r.terraflow, but try using
r.watershed.
Crashes with different errors. I'll post about that later. I was using terraflow because it claims to be for large data (not that I think this data is particularly large...)
r.watershed deals quite well with large datasets, too.See
https://grasswiki.osgeo.org/wiki/GRASS_GIS_Performance#Large_raster_data_processing
(r.watershed on a 90k x 100k pixels raster map)
r.mapcalc 'phi = bed*2' does produce a data file, full resolution or downsampled.
Oddly, "2*bed" works, but "2.5*bed" does not work.
That's rather unusual, indeed.
Does r.mapcalc only do integer math?
Definitely yes.
Not according to the manual.
Oh - so we should improve that:
https://grass.osgeo.org/grass70/manuals/r.mapcalc.html
Please let us know where to add extra explanations (what did you look
for and didn't find?)
"r.info bed" shows the range is -5028 to 3678.
Please note that if "bed" is an integer map, r.mapcalc will generate
integer from "bed * number". You would need to write
2.5 * bed
or
1.0 * bed * 2
to get floating point results.
Should we expand this section?
https://grass.osgeo.org/grass70/manuals/r.mapcalc.html#floating-point-values-in-the-expression
Markus
Hi Ken,
On Sun, Jan 3, 2016 at 6:06 PM, Markus Neteler <neteler@osgeo.org> wrote:
> r.mapcalc 'phi = bed*2' does produce a data file, full resolution or
downsampled.
>
> Oddly, "2*bed" works, but "2.5*bed" does not work.
That's rather unusual, indeed.
these expressions work for me:
r.mapcalc "test = 21"
r.mapcalc "test_a = 2*test"
r.mapcalc "test_b = 2.5*test"
r.mapcalc "test_c = test*2"
r.mapcalc "test_d = test*2.5"
r.info test_a | grep "Data Type"
| Data Type: CELL
r.info test_b | grep "Data Type"
| Data Type: DCELL
r.info test_c | grep "Data Type"
| Data Type: CELL
r.info test_d | grep "Data Type"
| Data Type: DCELL
CELL is for integer and DCELL for double (floating point). Spaces around
operators are best practice but nor required and I did not include them to
make it similar to you example.
> Does r.mapcalc only do integer math?
Definitely yes.
I'm not sure what Markus is referring ti but in general, r.mapcalc works
with integers and floating point numbers (floats and doubles) as requested.
> Not according to the manual.
Oh - so we should improve that:
https://grass.osgeo.org/grass70/manuals/r.mapcalc.html
Please let us know where to add extra explanations (what did you look
for and didn't find?)
In any case this is true, if you are/were confused about the functionality,
please suggest a change to the manual for benefit of others and yours as
well when you get back to it after some time.
Vaclav
On Jan 4, 2016 12:06 AM, “Markus Neteler” <neteler@osgeo.org> wrote:
On Sun, Jan 3, 2016 at 11:55 PM, Ken Mankoff <mankoff@gmail.com> wrote:
…
Does r.mapcalc only do integer math?
Definitely yes.
Ops, late night typo.
Correction: r.mapcalc does both integer and floating point math.
Markus
Hi,
I've figured out the issue, which is that I was looking at files on disk for some reason. I assume I shouldn't actually be in the /path/to/PERMANENT/ folder and use the x.info commands instead. Anyway, phi in the cell folder is 0 because it is a float, and the file does exist in the fcell folder.
On 2016-01-03 at 20:52, Vaclav Petras <wenzeslaus@gmail.com> wrote:
On Sun, Jan 3, 2016 at 6:06 PM, Markus Neteler <neteler@osgeo.org>
> Does r.mapcalc only do integer math?
> Not according to the manual.
Oh - so we should improve that:
https://grass.osgeo.org/grass70/manuals/r.mapcalc.html
Please let us know where to add extra explanations (what did you
look for and didn't find?)
The manual is clear. My double-negative sentence above was not clear.
Thanks for the helpful replies
-k.
P.S. I'm using grass in Emacs Org-Mode. I find it a nice setup. There is no feedback (that is an Org limitation), but it is really nice to be able to write reproducible documents and interleave text/comments/manuscript with the GIS commands.
On Mon, Jan 4, 2016 at 8:51 AM, Ken Mankoff <mankoff@gmail.com> wrote:
I've figured out the issue, which is that I was looking at files on disk
for some reason. I assume I shouldn't actually be in the
/path/to/PERMANENT/ folder and use the x.info commands instead. Anyway,
phi in the cell folder is 0 because it is a float, and the file does exist
in the fcell folder.
Yes, the file in cell dir is there for each raster. And yes, r.info is the
right thing. GRASS Database, Location and Mapset are database, so special
tools like r.info are required.
P.S. I'm using grass in Emacs Org-Mode. I find it a nice setup. There is
no feedback (that is an Org limitation), but it is really nice to be able
to write reproducible documents and interleave text/comments/manuscript
with the GIS commands.
Interesting. Do you plan any blog post about it?
On 2016-01-04 at 10:18, Vaclav Petras <wenzeslaus@gmail.com> wrote:
Interesting. Do you plan any blog post about it?
I will, but not until after the paper is published (which will be a while...). I plan to include the Org file as supplemental material so readers can fully reproduce it.
-k.