[GRASS-user] A question before I embark on a programming exercise

GRASS GRASS user
1

All:

I need to generate an ascii text file from a flow direction grid that consists of (among a couple other things that don’t really matter at this point) for each pixel:

(1) a unique integer identifier (1 – N) for the pixel

(2) the integer identifier of the downstream pixel (assuming there is ONLY one)

(3) the x,y location of the pixel (presumably, the lower left corner of the pixel)

has anyone already done something like this?

It’s needed (along with some header information) as an input file for a gridded distributed hydrologic model, identifying the flow connectivity from pixel to pixel for streamflow routing purposes. I have already formulated conceptually how to do this, but if someone has already done such a thing, why reinvent the wheel??

Any help is appreciated.

Regards,

Tom

2

Thomas Adams <tea3rd@gmail.com> writes:

All:

I need to generate an ascii text file from a flow direction grid that
consists of (among a couple other things that don't really matter at this
point) for each pixel:

(1) a unique integer identifier (1 -- N) for the pixel
(2) the integer identifier of the downstream pixel (assuming there is ONLY
one)
(3) the x,y location of the pixel (presumably, the lower left corner of the
pixel)

has anyone already done something like this?

Well - I have done something similar, i.e. simulated water dispersal, in
a package which is at [1]. The actual function is at [1]. After reading
your mail again, it seems to be different, but it still might give you
some ideas.

Cheers,

Rainer

It's needed (along with some header information) as an input file for a
gridded distributed hydrologic model, identifying the flow connectivity
from pixel to pixel for streamflow routing purposes. I have already
formulated conceptually how to do this, but if someone has already done
such a thing, why reinvent the wheel??

Any help is appreciated.

Regards,
Tom
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Footnotes:
[1] https://github.com/rkrug/seedDisp

[1] https://github.com/rkrug/seedDisp/blob/master/R/waterDispGRASS.R

--
Rainer M. Krug
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982

3

Rainer,

Thank you very much. I just downloaded everything via git. I greatly appreciate your help.

Best,

Tom

···

On Wed, Nov 12, 2014 at 1:28 AM, Rainer M Krug <Rainer@krugs.de> wrote:

Thomas Adams <tea3rd@gmail.com> writes:

All:

I need to generate an ascii text file from a flow direction grid that
consists of (among a couple other things that don’t really matter at this
point) for each pixel:

(1) a unique integer identifier (1 – N) for the pixel
(2) the integer identifier of the downstream pixel (assuming there is ONLY
one)
(3) the x,y location of the pixel (presumably, the lower left corner of the
pixel)

has anyone already done something like this?

Well - I have done something similar, i.e. simulated water dispersal, in
a package which is at [1]. The actual function is at [1]. After reading
your mail again, it seems to be different, but it still might give you
some ideas.

Cheers,

Rainer

It’s needed (along with some header information) as an input file for a
gridded distributed hydrologic model, identifying the flow connectivity
from pixel to pixel for streamflow routing purposes. I have already
formulated conceptually how to do this, but if someone has already done
such a thing, why reinvent the wheel??

Any help is appreciated.

Regards,
Tom

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Footnotes:
[1] https://github.com/rkrug/seedDisp

[1] https://github.com/rkrug/seedDisp/blob/master/R/waterDispGRASS.R

–
Rainer M. Krug
email: Rainerkrugsde
PGP: 0x0F52F982

4

On 12/11/14 05:10, Thomas Adams wrote:

I need to generate an ascii text file from a flow direction grid that
consists of (among a couple other things that don't really matter at
this point) for each pixel:

Rainer has already given you a whole program, but here's a decomposition in simple (untested) r.mapcalc steps:

(1) a unique integer identifier (1 -- N) for the pixel

r.mapcalc "id = row()*N + col()" # where N = nb cols * 10

(you might have to watch out for integer overflow, though)

(2) the integer identifier of the downstream pixel (assuming there is
ONLY one)

Calculate flow direction in SFD mode and use something like this:

r.mapcalc down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3, id[-1,1] .... etc

(3) the x,y location of the pixel (presumably, the lower left corner of
the pixel)

center of pixel:
r.mapcalc x = x()
r.mapcalc y = y()

if you want the lower left corner:

r.mapcalc x = x() - ewres()/2
r.mapcalc y = y() - nsres()/2

Then export everything with

r.stats -1 in=id,x,y,down_id

Moritz

5

Awesome, Moritz!

Thank you so much…

Tom

···

On Wed, Nov 12, 2014 at 7:30 AM, Moritz Lennert <mlennert@club.worldonline.be> wrote:

On 12/11/14 05:10, Thomas Adams wrote:

I need to generate an ascii text file from a flow direction grid that
consists of (among a couple other things that don’t really matter at
this point) for each pixel:

Rainer has already given you a whole program, but here’s a decomposition in simple (untested) r.mapcalc steps:

(1) a unique integer identifier (1 – N) for the pixel

r.mapcalc “id = row()*N + col()” # where N = nb cols * 10

(you might have to watch out for integer overflow, though)

(2) the integer identifier of the downstream pixel (assuming there is
ONLY one)

Calculate flow direction in SFD mode and use something like this:

r.mapcalc down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3, id[-1,1] … etc

(3) the x,y location of the pixel (presumably, the lower left corner of
the pixel)

center of pixel:
r.mapcalc x = x()
r.mapcalc y = y()

if you want the lower left corner:

r.mapcalc x = x() - ewres()/2
r.mapcalc y = y() - nsres()/2

Then export everything with

r.stats -1 in=id,x,y,down_id

Moritz

6

Hi Tom,

what comes into my mind is following approach (not tested):

  1. calculate flow direction using r.watershed for each pixel of your map
  2. transforming raster to vector points (r.to.vect), adding some columns for your results (X, Y, ID). This step should provide you already with unique IDs for your cells/points.
  3. With v.to.db you can update the X and Y coordinate of your cell center and together with the raster resolution you can get the lower left corner as desired.
  4. By knowing the flow direction for each point you should be able to calculate the X and Y coordinate of the downstream cell (adding or substracting your cell dimensions from the cell center, correct for diagonal flow). When you know X and Y of your downstream cell then you should be able to also get the ID of the downstream cell (most probably a SQL query I guess)

Maybe there are easier ways…

hope it works,
Johannes

···

On Wed, Nov 12, 2014 at 5:10 AM, Thomas Adams <tea3rd@gmail.com> wrote:

All:

I need to generate an ascii text file from a flow direction grid that consists of (among a couple other things that don’t really matter at this point) for each pixel:

(1) a unique integer identifier (1 – N) for the pixel

(2) the integer identifier of the downstream pixel (assuming there is ONLY one)

(3) the x,y location of the pixel (presumably, the lower left corner of the pixel)

has anyone already done something like this?

It’s needed (along with some header information) as an input file for a gridded distributed hydrologic model, identifying the flow connectivity from pixel to pixel for streamflow routing purposes. I have already formulated conceptually how to do this, but if someone has already done such a thing, why reinvent the wheel??

Any help is appreciated.

Regards,

Tom

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

7

Johannes,

Thank you very much for the useful suggestions. I greatly appreciate your comments! I’ll post my results once I work through things.

Best,

Tom

···

On Wed, Nov 12, 2014 at 7:42 AM, Johannes Radinger <johannesradinger@gmail.com> wrote:

Hi Tom,

what comes into my mind is following approach (not tested):

  1. calculate flow direction using r.watershed for each pixel of your map
  2. transforming raster to vector points (r.to.vect), adding some columns for your results (X, Y, ID). This step should provide you already with unique IDs for your cells/points.
  3. With v.to.db you can update the X and Y coordinate of your cell center and together with the raster resolution you can get the lower left corner as desired.
  4. By knowing the flow direction for each point you should be able to calculate the X and Y coordinate of the downstream cell (adding or substracting your cell dimensions from the cell center, correct for diagonal flow). When you know X and Y of your downstream cell then you should be able to also get the ID of the downstream cell (most probably a SQL query I guess)

Maybe there are easier ways…

hope it works,
Johannes

–

On Wed, Nov 12, 2014 at 5:10 AM, Thomas Adams <tea3rd@gmail.com> wrote:

All:

I need to generate an ascii text file from a flow direction grid that consists of (among a couple other things that don’t really matter at this point) for each pixel:

(1) a unique integer identifier (1 – N) for the pixel

(2) the integer identifier of the downstream pixel (assuming there is ONLY one)

(3) the x,y location of the pixel (presumably, the lower left corner of the pixel)

has anyone already done something like this?

It’s needed (along with some header information) as an input file for a gridded distributed hydrologic model, identifying the flow connectivity from pixel to pixel for streamflow routing purposes. I have already formulated conceptually how to do this, but if someone has already done such a thing, why reinvent the wheel??

Any help is appreciated.

Regards,

Tom

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Thomas E Adams, III
718 McBurney Drive
Lebanon, OH 45036

1 (513) 739-9512 (cell)

8

Moritz,

I need to be sure I understand this line:

down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3, id[-1,1] … etc

dir=1 is what direction?; dir=2 is what direction?, etc. Is dir=1 directly to the ‘East’, dir=3 directly ‘south’, dir=5 directly ‘west’ and dir=7 directly ‘north’

Also, does the indexing [0,1] refer to the [x,y], i.e., [column, row] so that [0,1] refers to the same column, but one row below? Which would imply from your line that dir=2 is directly ‘south’??

I apologize for being dense about this, I just need to be sure — sorry…

I did get the first part to work using r.mapcalc to get the pixel indexing: id = (row() - 1)*Ncolumns + col()

Regards,
Tom

···

On Wed, Nov 12, 2014 at 7:30 AM, Moritz Lennert <mlennert@club.worldonline.be> wrote:

On 12/11/14 05:10, Thomas Adams wrote:

I need to generate an ascii text file from a flow direction grid that
consists of (among a couple other things that don’t really matter at
this point) for each pixel:

Rainer has already given you a whole program, but here’s a decomposition in simple (untested) r.mapcalc steps:

(1) a unique integer identifier (1 – N) for the pixel

r.mapcalc “id = row()*N + col()” # where N = nb cols * 10

(you might have to watch out for integer overflow, though)

(2) the integer identifier of the downstream pixel (assuming there is
ONLY one)

Calculate flow direction in SFD mode and use something like this:

r.mapcalc down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3, id[-1,1] … etc

(3) the x,y location of the pixel (presumably, the lower left corner of
the pixel)

center of pixel:
r.mapcalc x = x()
r.mapcalc y = y()

if you want the lower left corner:

r.mapcalc x = x() - ewres()/2
r.mapcalc y = y() - nsres()/2

Then export everything with

r.stats -1 in=id,x,y,down_id

Moritz

9

On 13/11/14 00:25, Thomas Adams wrote:

Moritz,

I need to be sure I understand this line:

down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3, id[-1,1] .... etc

dir=1 is what direction?; dir=2 is what direction?, etc. Is dir=1
directly to the 'East', dir=3 directly 'south', dir=5 directly 'west'
and dir=7 directly 'north'

In the example, I was assuming (without making it explicit, sorry) a dir coming from r.watershed ('drainage' parameter) and as you can see in the man page:

"Provides the "aspect" for each cell measured CCW from East. Multiplying positive values by 45 will give the direction in degrees that the surface runoff will travel from that cell. The value 0 (zero) indicates that the cell is a depression area (defined by the depression input map). Negative values indicate that surface runoff is leaving the boundaries of the current geographic region. The absolute value of these negative cells indicates the direction of flow."

I.e. 2 = due north, 4 = east, 6 = south, 8 = east.

But obviously this has to be adapted to whatever tool you use for calculating your drainage direction.

Also, does the indexing [0,1] refer to the [x,y], i.e., [column, row] so
that [0,1] refers to the same column, but one row below? Which would
imply from your line that dir=2 is directly 'south'??

From the r.mapcalc man page:

"In r.mapcalc, maps may be followed by a neighborhood modifier that specifies a relative offset from the current cell being evaluated. The format is map[r,c], where r is the row offset and c is the column offset. For example, map[1,2] refers to the cell one row below and two columns to the right of the current cell, map[-2,-1] refers to the cell two rows above and one column to the left of the current cell, and map[0,1] refers to the cell one column to the right of the current cell."

I.e. the inverse of what you interpreted: [0,1] = same row, one to the east, i.e. drainage direction due east.

And I realise that I inversed it myself in the second element of the example. It should be:

down_id = if(dir=1, id[1,1], if(dir=2, id[1,0], if(dir=3, id[-1,1] .... etc

i.e. dir=2 implies due North and so one row up, same column.

I apologize for being dense about this, I just need to be sure — sorry…

It can be a bit difficult to wrap your head around, especially when those trying to help you induce you into error :wink:

Moritz

10

Moritz Lennert wrote:

r.mapcalc down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3,
id[-1,1] .... etc

r.mapcalc follows C expression syntax; "="is assignment, "==" is the
equality test.

--
Glynn Clements <glynn@gclements.plus.com>

11

Glynn Clements <glynn@gclements.plus.com> writes:

Moritz Lennert wrote:

r.mapcalc down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3,
id[-1,1] .... etc

r.mapcalc follows C expression syntax; "="is assignment, "==" is the
equality test.

In this case shouldn't it be

,----
| r.mapcalc down_id = if(dir==1, id[1,1], if(dir==2, id[0,1], if(dir==3, id[-1,1] .... etc
`----

with "==" after the dir?

Rainer

--
Rainer M. Krug
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982

12

On 14/11/14 11:21, Rainer M Krug wrote:

Glynn Clements <glynn@gclements.plus.com> writes:

Moritz Lennert wrote:

r.mapcalc down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3,
id[-1,1] .... etc

r.mapcalc follows C expression syntax; "="is assignment, "==" is the
equality test.

In this case shouldn't it be

,----
| r.mapcalc down_id = if(dir==1, id[1,1], if(dir==2, id[0,1], if(dir==3, id[-1,1] .... etc
`----

with "==" after the dir?

Yes, that's what Glynn is saying. Sorry, I wrote this example too quickly.

Moritz

13

Glynn & Rainer,

YES, absolutely – I caught that and made it ‘==’ instead of ‘=’. When typing quickly, I sometimes make that mistake and then have to correct myself. But fundamentally, Moritz’ approach worked very well and produced what I needed. I have a more difficult problem now, something I think I can only resolve by using higher resolution data. But, that is OK, too, because I should get better modeling results in the end. I’m trying to create my own approach to generating the file, that when used by the model, gives essentially the same results.

The folks who wrote the model and software to generate this runoff/flow connectivity file between pixels, found that with coarser pixel resolutions, the connection between pixels was badly distorted, in physically unreasonable ways. I don’t have their original code, so I can not use it, plus, I really want to model at much finer spatial resolutions anyway, where the problem should go vanish.

So, I just need to use higher resolution DEM data, which I have, rather than trying to do things at the default 4km resolution.

One question I have: it would be very nice to be able to show how the pixels are connected using vectors, essentially the drainage network. Any thoughts how to do this?

Regards,

Tom

···

On Fri, Nov 14, 2014 at 3:21 AM, Rainer M Krug <Rainer@krugs.de> wrote:

Glynn Clements <glynn@gclements.plus.com> writes:

Moritz Lennert wrote:

r.mapcalc down_id = if(dir=1, id[1,1], if(dir=2, id[0,1], if(dir=3,
id[-1,1] … etc

r.mapcalc follows C expression syntax; "="is assignment, “==” is the
equality test.

In this case shouldn’t it be

,----
| r.mapcalc down_id = if(dir==1, id[1,1], if(dir==2, id[0,1], if(dir==3, id[-1,1] … etc
`----

with “==” after the dir?

Rainer

–
Rainer M. Krug
email: Rainerkrugsde
PGP: 0x0F52F982

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user