Hi list,
I'm now trying to use r.watershed instead of r.terraflow. It works with my surface DEM, but I want to use a different 'cost' function as my surface instead of the actual land surface. I generate this cost function with:
r.mapcalc 'phi = (1000 * 9.8 * bed + 917 * 9.8 * (surf - bed))'
Even before r.watershed (which works on bed and surf but not on phi), I notice phi is different in two ways:
1) phi has 'categories'
r.info bed | grep Cat
r.info surf | grep Cat
r.info phi | grep Cat
| Type of Map: raster Number of Categories: 0 |
| Type of Map: raster Number of Categories: 0 |
| Type of Map: raster Number of Categories: 255 |
2) when I display (d.erase; d.rast bed; d.legend map=bed) bed or surf, the colorbar is green->blue->red (top-mid-bottom) and numbers are low->high (top->bottom). When I display phi, the order is reversed.
When I run r.watershed, it appears to work, but the results are garbage. flow directions are all E or W (not 8-directions as the same command generates w/ surf or bed), etc. But I assume these r.watershed issues are secondary.
Can someone help clarify what is going on? I'm happy to provide a full script w/ wget to access the data if it helps.
Thanks,
-k.
Ken,
On Tue, Jan 5, 2016 at 3:28 AM, Ken Mankoff <mankoff@gmail.com> wrote:
...
Can someone help clarify what is going on? I'm happy to provide a full script w/ wget to access the data if it helps.
This would be very useful.
Markus
Hi Markus,
Thanks for the reply.
On 2016-01-06 at 16:25, Markus Neteler <neteler@osgeo.org> wrote:
On Tue, Jan 5, 2016 at 3:28 AM, Ken Mankoff <mankoff@gmail.com> wrote:
...
Can someone help clarify what is going on? I'm happy to provide a
full script w/ wget to access the data if it helps.
This would be very useful.
Before I supply a full script with external data, I wonder if I can make synthetic data to debug/recreate this more easily. I see that I can do "r.mapcalc 'x = 42'" to create a new vector raster with 1 point. Is there some easy method to create a 2D raster within grass? Or should I make a simple ASCII CSV file and read that in? If I can recreate the behavior with a small hand-made data set I think it would help me figure out what is wrong with the data.
Thanks,
-k.
Hi Ken,
…
Before I supply a full script with external data, I wonder if I can make synthetic data to debug/recreate this more easily. I see that I can do “r.mapcalc ‘x = 42’” to create a new vector raster with 1 point. Is there some easy method to create a 2D raster within grass? Or should I make a simple ASCII CSV file and read that in? If I can recreate the behavior with a small hand-made data set I think it would help me figure out what is wrong with the data.
r.mapcalc uses the region definition to create the new maps. So, if your region is one cell (you can check with g.region -p) you’ll have a one cell raster as output. Therefore, if you want a 2D raster, just set a bigger region, something like g.region cols=2 rows=2 (check the manual for more options to set the region), and then the r.mapcalc operation you need.
Hope it helps,
Vero
On Wed, Jan 6, 2016 at 9:57 PM, Ken Mankoff <mankoff@gmail.com> wrote:
Is there some easy method to create a 2D raster within grass?
I’m not sure what you need but you can also create a nice surface using r.surf.fractal.
https://grass.osgeo.org/grass70/manuals/r.surf.fractal.html
Hi,
So while preparing a minimal example the bug got resolved. I think it had to do with stale variables (or a stale mask) that remained in the session. It seems the cause wasn't in the software but something between the keyboard and monitor outside of the laptop...
Thanks for the help so far,
-k.
On 2016-01-06 at 22:46, Veronica Andreo <veroandreo@gmail.com> wrote:
Hi Ken,
> ...
Before I supply a full script with external data, I wonder if I can
make synthetic data to debug/recreate this more easily. I see that I
can do "r.mapcalc 'x = 42'" to create a new vector raster with 1
point. Is there some easy method to create a 2D raster within grass?
Or should I make a simple ASCII CSV file and read that in? If I can
recreate the behavior with a small hand-made data set I think it would
help me figure out what is wrong with the data.
r.mapcalc uses the region definition to create the new maps. So, if
your region is one cell (you can check with g.region -p) you'll have a
one cell raster as output. Therefore, if you want a 2D raster, just
set a bigger region, something like g.region cols=2 rows=2 (check the
manual for more options to set the region), and then the r.mapcalc
operation you need.
Hope it helps,
Vero