[GRASS-dev] report on raster point creation

[cmbarton]

Hamish and Helena,

We have finally had an initial production test of methods to create a raster point using xy coordinates. If you remember, I’ve periodically whined a bit about the lack of a way to create and change a raster cell based on its xy coordinates instead of its cat value.

The 2 suggestions I got to work around this were to…

1. use r.in.poly to create the raster map; make an input file such that for each point, the input is set to line (L) with 2 points that are identical.
2. use v.in.ascii and then use v.to.rast to create the raster map.

We are coupling an agent-based simulation to GRASS, such that the simulated land use activities will affect land cover, which in turn will affect a landscape evolution routine.

Optimizing the speed of this is critical, and so we prefer to reduce the number of steps needed to accomplish this. This means that the r.in.poly approach is the one we prefer. However, in tests conducted today, we ran into an unexpected problem.

If the coordinate of the point to be created happens to fall on a cell boundary, we get a 2-cell ‘line’ instead of a point even though the beginning and end points are identical. Apparently, r.in.poly looks for a cell center closest to end point one and then searches for a cell center for end point 2 that is different from the first cell, resulting in a 2 cell line. This does not happen if the point does not fall on a cell boundary. Attached is the raster file created from the following r.in.poly input. The first input creates the 2 yellow cells, the other 2 create the red and green single cells.

L
592250.0 4924400.0
592250.0 4924400.0
= 40 household0

L
592362.0 4924531.0
592362.0 4924531.0
= 41 household1

L
592142.0 4924277.0
592142.0 4924277.0
= 43 household3

This is a potential problem, since an agent might not know that the the coordinates it chooses for a field to farm or graze like on a cell boundary. It is OK if r.in.poly substitutes the nearest grid cell center to make the raster point for the xy coordinates that the agent chooses, but not OK if it picks 2 different cells instead of one.

The simplest solution is to allow r.in.poly to make points as well as lines and areas. I’m not sure why it doesn’t allow points in the first place. It seems like an oversight. Is this difficult to implement?

---

Also, there is still no one-step way to change a cell value based on its xy coordinates. You can change it based on its position relative to another cell, and based on its cat, but not its coordinates. Currently, the only way to do this is to create some cells (e.g., using r.in.poly) and then patch them back into the original map, overwriting it. This is not a problem for our own current modeling efforts, given the way we are doing it (making an impact map of agent landuse activities and using r.mapcalc to alter a land cover map accordingly), but it might help us better optimize the running of these models in the future.

However, if there is interest in doing agent and cellular modeling within a GRASS platform, it will become necessary to be able to do this to simulate dynamics. Patching may simply not be practical for this kind of work. As Python becomes more used as a scripting environment, this is an attractive kind of application to build in Python, given its object oriented structure, and couple with GRASS. Looking ahead, I can see such modeling as having being increasingly important to researchers doing research in the social and natural sciences. Coupled with an OO language like Python, GRASS could be a real leading application in this area. So another whine from me about this.

Michael

---

Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

(attachments)

[H_B]

Michael Barton wrote:

We have finally had an initial production test of methods to create a
raster point using xy coordinates. If you remember, I've periodically
whined a bit about the lack of a way to create and change a raster
cell based on its xy coordinates instead of its cat value.

The 2 suggestions I got to work around this were to...

1) use r.in.poly to create the raster map; make an input file such
that for each point, the input is set to line (L) with 2 points that
are identical.

To make that method easier I have just added a "P" instruction to
r.in.poly in 6.3cvs. That uses G_plot_point(), which draws a 0 length
line using the same method, so probably it doesn't fix your 2 cell if
on bound problem (?? try it and see).

2) use v.in.ascii and then use v.to.rast to create the raster map.

Note the v.to.rast code is ~95% a clone of the r.in.poly code, so results
(and problems) should be mostly the same. Do you see the same multiple
cell trouble if the point falls on a bound using this method?

It probably needs to be noted in the v.to.rast help page if it does, at
minimum.

3) use r.in.xyz (I expect this will be the fastest method)

4) r.mapcalc outmap='if( (abs(x() - $x)) < ewres()) \
      && (abs(y() - $y) < nsres(), $val, null() )'
[or something like that, for a single point; yuk]

5) do the gridding yourself and create a r.in.ascii file.

all these would need to be followed by a call to r.patch.
(see d.rast.edit.tcl method)

We are coupling an agent-based simulation to GRASS, such that the
simulated land use activities will affect land cover, which in turn
will affect a landscape evolution routine.

are you updating a sequential string of cells, blocks of cells, or
random cells at each update?

Optimizing the speed of this is critical, and so we prefer to reduce
the number of steps needed to accomplish this. This means that the
r.in.poly approach is the one we prefer. However, in tests conducted
today, we ran into an unexpected problem.

If the coordinate of the point to be created happens to fall on a cell
boundary, we get a 2-cell 'line' instead of a point even though the
beginning and end points are identical. Apparently, r.in.poly looks
for a cell center closest to end point one and then searches for a
cell center for end point 2 that is different from the first cell,
resulting in a 2 cell line. This does not happen if the point does not
fall on a cell boundary. Attached is the raster file created from the
following r.in.poly input. The first input creates the 2 yellow cells,
the other 2 create the red and green single cells.

r.in.poly/raster.c defines the "cont" continue line fn used by
G_setup_plot() (and thus all following G_plot_*() commands) with
G_bresenham_line(). That in turn uses the "dot" fn to draw the points;
which "dot" is used depends on the result of the getformat() scan of
the cat values (& I'm about foggy about what that's all about).

for info on the Bresenham line rasterization method see:
  http://en.wikipedia.org/wiki/Bresenham’s_line_algorithm

This is a potential problem, since an agent might not know that the
the coordinates it chooses for a field to farm or graze like on a cell
boundary. It is OK if r.in.poly substitutes the nearest grid cell
center to make the raster point for the xy coordinates that the agent
chooses, but not OK if it picks 2 different cells instead of one.

Try r.in.xyz, it shouldn't suffer from this: on one side of the cell it
tests for <= and on the other side <. So a point on a bound always fall
in the <= cell. The exception is the outer E,W raster boundaries which
will make values fall "into" the raster if on the bound.

r.in.xyz doesn't let you add cat labels, but you can do that manually by
creating a cats/ file. See spearfish60/PERMANENT/cats/landuse for an
example, it's very simple.

The simplest solution is to allow r.in.poly to make points as well as
lines and areas. I'm not sure why it doesn't allow points in the first
place. It seems like an oversight. Is this difficult to implement?

done, but I don't think it will help any due to the way G_plot_point()
works.

* * * *

Also, there is still no one-step way to change a cell value based on
its xy coordinates.

So you want a non-interactive d.rast.edit, like v.edit?

A couple of ways to do that, 1) get input x1,y1,cat1 ... from stdin,
qsort() based on y then x, read each line and if there's an updated
cell in it replace it before writing the row back out. (as r.patch)
2) calculate the array column,row which contains each x,y and update
each in the order they arrive (as r.in.xyz). Probably the fastest way is
a combination of the two. Sort by y, and for each row update the column
values in the order they arrive, then write out the finished line and
move on to the next row.

However, if there is interest in doing agent and cellular modeling
within a GRASS platform, it will become necessary to be able to do
this to simulate dynamics. Patching may simply not be practical for
this kind of work. As Python becomes more used as a scripting
environment, this is an attractive kind of application to build in
Python, given its object oriented structure, and couple with GRASS.

The module probably isn't so hard to write, in C. r.example should give
a nice start.

It would be interesting to rewrite r.example.c using the SWIG Python
interface, and see if that gains any simplicity over C, or if it just
adds complexity without simplicity -- as you still need to go through
all the same chores.

Hamish

[cmbarton]

Thanks very much. I'll need to digest all this useful information later
today, but it is helpful. I'll check r.in.xyz and your update to r.in.poly
and get back to you.

Michael

On 6/14/07 2:45 AM, "Hamish" <hamish_nospam@yahoo.com> wrote:

Michael Barton wrote:

We have finally had an initial production test of methods to create a
raster point using xy coordinates. If you remember, I've periodically
whined a bit about the lack of a way to create and change a raster
cell based on its xy coordinates instead of its cat value.

The 2 suggestions I got to work around this were to...

1) use r.in.poly to create the raster map; make an input file such
that for each point, the input is set to line (L) with 2 points that
are identical.

To make that method easier I have just added a "P" instruction to
r.in.poly in 6.3cvs. That uses G_plot_point(), which draws a 0 length
line using the same method, so probably it doesn't fix your 2 cell if
on bound problem (?? try it and see).

2) use v.in.ascii and then use v.to.rast to create the raster map.

Note the v.to.rast code is ~95% a clone of the r.in.poly code, so results
(and problems) should be mostly the same. Do you see the same multiple
cell trouble if the point falls on a bound using this method?

It probably needs to be noted in the v.to.rast help page if it does, at
minimum.

3) use r.in.xyz (I expect this will be the fastest method)

4) r.mapcalc outmap='if( (abs(x() - $x)) < ewres()) \
      && (abs(y() - $y) < nsres(), $val, null() )'
[or something like that, for a single point; yuk]

5) do the gridding yourself and create a r.in.ascii file.

all these would need to be followed by a call to r.patch.
(see d.rast.edit.tcl method)

We are coupling an agent-based simulation to GRASS, such that the
simulated land use activities will affect land cover, which in turn
will affect a landscape evolution routine.

are you updating a sequential string of cells, blocks of cells, or
random cells at each update?

Optimizing the speed of this is critical, and so we prefer to reduce
the number of steps needed to accomplish this. This means that the
r.in.poly approach is the one we prefer. However, in tests conducted
today, we ran into an unexpected problem.

If the coordinate of the point to be created happens to fall on a cell
boundary, we get a 2-cell 'line' instead of a point even though the
beginning and end points are identical. Apparently, r.in.poly looks
for a cell center closest to end point one and then searches for a
cell center for end point 2 that is different from the first cell,
resulting in a 2 cell line. This does not happen if the point does not
fall on a cell boundary. Attached is the raster file created from the
following r.in.poly input. The first input creates the 2 yellow cells,
the other 2 create the red and green single cells.

r.in.poly/raster.c defines the "cont" continue line fn used by
G_setup_plot() (and thus all following G_plot_*() commands) with
G_bresenham_line(). That in turn uses the "dot" fn to draw the points;
which "dot" is used depends on the result of the getformat() scan of
the cat values (& I'm about foggy about what that's all about).

for info on the Bresenham line rasterization method see:
  Bresenham's line algorithm - Wikipedia

This is a potential problem, since an agent might not know that the
the coordinates it chooses for a field to farm or graze like on a cell
boundary. It is OK if r.in.poly substitutes the nearest grid cell
center to make the raster point for the xy coordinates that the agent
chooses, but not OK if it picks 2 different cells instead of one.

Try r.in.xyz, it shouldn't suffer from this: on one side of the cell it
tests for <= and on the other side <. So a point on a bound always fall
in the <= cell. The exception is the outer E,W raster boundaries which
will make values fall "into" the raster if on the bound.

r.in.xyz doesn't let you add cat labels, but you can do that manually by
creating a cats/ file. See spearfish60/PERMANENT/cats/landuse for an
example, it's very simple.

The simplest solution is to allow r.in.poly to make points as well as
lines and areas. I'm not sure why it doesn't allow points in the first
place. It seems like an oversight. Is this difficult to implement?

done, but I don't think it will help any due to the way G_plot_point()
works.

* * * *

Also, there is still no one-step way to change a cell value based on
its xy coordinates.

So you want a non-interactive d.rast.edit, like v.edit?

A couple of ways to do that, 1) get input x1,y1,cat1 ... from stdin,
qsort() based on y then x, read each line and if there's an updated
cell in it replace it before writing the row back out. (as r.patch)
2) calculate the array column,row which contains each x,y and update
each in the order they arrive (as r.in.xyz). Probably the fastest way is
a combination of the two. Sort by y, and for each row update the column
values in the order they arrive, then write out the finished line and
move on to the next row.

However, if there is interest in doing agent and cellular modeling
within a GRASS platform, it will become necessary to be able to do
this to simulate dynamics. Patching may simply not be practical for
this kind of work. As Python becomes more used as a scripting
environment, this is an attractive kind of application to build in
Python, given its object oriented structure, and couple with GRASS.

The module probably isn't so hard to write, in C. r.example should give
a nice start.

It would be interesting to rewrite r.example.c using the SWIG Python
interface, and see if that gains any simplicity over C, or if it just
adds complexity without simplicity -- as you still need to go through
all the same chores.

Hamish

__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

[David_Finlayson2]

Michael -

I often couple GRASS to numerical models or complex grid-based computations using the following method:

Step 1: Export the environmental grids needed for the model using r.out.arc (an easy format to parse in Python) or a simple binary format.
Step 2: Read the data into a Python class that represents the grid as a 2D array (numeric is great for this). The class also has a few helper methods such as: read_grid(), write_grid(), get(x, y), set(x,y), etc.
Step 3: Run the model
Step 4: Dump the python grid back to the ASCII Grid format (or raw binary)
Step 5: Read the data back into GRASS using r.in.arc (etc)

Most of the work in the model occurs in Steps 2-4 where there is no overhead on the model due to carrying around the GIS baggage. You could use python, R, Matlab, whatever.

As an unsolicited aside, C is actually a pretty good language for grid-based models when one uses a numerical library such as the GNU Scientific Library to handle the matrix (grid) stuff. Two bonuses to using C are the increase in performance of the model and by the end of the project, you have learned enough C to read and start using the GRASS source code directly. Not only that, but there are lots of GIS-related libraries ( i.e., GDAL, Proj4, Shapelib) that are easy to use from C once you learn the basics.

I resisted C for years because Python was good enough and C was “hard”. Today I think that was being penny-wise and pound foolish. There are many projects in GIS where speed really matters. Whether you are filtering a very large Lidar data set, or running a very small model millions of times (such as in a Monte Carlo simulation), speed really matters in these cases. The huge library of existing C code is real boon, too. I find that for small projects (< 10k lines of code) C is a pretty good language; its warts don’t show up until project get really large. And like I said, it’s nice to be able to customize GRASS source when you need to.

Cheers,

David

On 6/14/07, Michael Barton <michael.barton@asu.edu> wrote:

Thanks very much. I’ll need to digest all this useful information later
today, but it is helpful. I’ll check r.in.xyz and your update to r.in.poly
and get back to you.

Michael

On 6/14/07 2:45 AM, “Hamish” <hamish_nospam@yahoo.com> wrote:

Michael Barton wrote:

We have finally had an initial production test of methods to create a
raster point using xy coordinates. If you remember, I’ve periodically
whined a bit about the lack of a way to create and change a raster
cell based on its xy coordinates instead of its cat value.

The 2 suggestions I got to work around this were to…

1. use r.in.poly to create the raster map; make an input file such
   that for each point, the input is set to line (L) with 2 points that
   are identical.

To make that method easier I have just added a “P” instruction to
r.in.poly in 6.3cvs. That uses G_plot_point(), which draws a 0 length
line using the same method, so probably it doesn’t fix your 2 cell if
on bound problem (?? try it and see).

2. use v.in.ascii and then use v.to.rast to create the raster map.

Note the v.to.rast code is ~95% a clone of the r.in.poly code, so results
(and problems) should be mostly the same. Do you see the same multiple
cell trouble if the point falls on a bound using this method?

It probably needs to be noted in the v.to.rast help page if it does, at
minimum.

3. use r.in.xyz (I expect this will be the fastest method)

4. r.mapcalc outmap=‘if( (abs(x() - $x)) < ewres())
   && (abs(y() - $y) < nsres(), $val, null() )’
   [or something like that, for a single point; yuk]

5. do the gridding yourself and create a r.in.ascii file.

all these would need to be followed by a call to r.patch.
(see d.rast.edit.tcl method)

We are coupling an agent-based simulation to GRASS, such that the
simulated land use activities will affect land cover, which in turn
will affect a landscape evolution routine.

are you updating a sequential string of cells, blocks of cells, or
random cells at each update?

Optimizing the speed of this is critical, and so we prefer to reduce
the number of steps needed to accomplish this. This means that the
r.in.poly approach is the one we prefer. However, in tests conducted
today, we ran into an unexpected problem.

If the coordinate of the point to be created happens to fall on a cell
boundary, we get a 2-cell ‘line’ instead of a point even though the
beginning and end points are identical. Apparently, r.in.poly looks
for a cell center closest to end point one and then searches for a
cell center for end point 2 that is different from the first cell,
resulting in a 2 cell line. This does not happen if the point does not
fall on a cell boundary. Attached is the raster file created from the
following r.in.poly input. The first input creates the 2 yellow cells,
the other 2 create the red and green single cells.

r.in.poly/raster.c defines the “cont” continue line fn used by
G_setup_plot() (and thus all following G_plot_*() commands) with
G_bresenham_line(). That in turn uses the “dot” fn to draw the points;
which “dot” is used depends on the result of the getformat() scan of
the cat values (& I’m about foggy about what that’s all about).

for info on the Bresenham line rasterization method see:
http://en.wikipedia.org/wiki/Bresenham’s_line_algorithm

This is a potential problem, since an agent might not know that the
the coordinates it chooses for a field to farm or graze like on a cell
boundary. It is OK if r.in.poly substitutes the nearest grid cell
center to make the raster point for the xy coordinates that the agent
chooses, but not OK if it picks 2 different cells instead of one.

Try r.in.xyz, it shouldn’t suffer from this: on one side of the cell it
tests for <= and on the other side <. So a point on a bound always fall
in the <= cell. The exception is the outer E,W raster boundaries which
will make values fall “into” the raster if on the bound.

r.in.xyz doesn’t let you add cat labels, but you can do that manually by
creating a cats/ file. See spearfish60/PERMANENT/cats/landuse for an
example, it’s very simple.

The simplest solution is to allow r.in.poly to make points as well as
lines and areas. I’m not sure why it doesn’t allow points in the first
place. It seems like an oversight. Is this difficult to implement?

done, but I don’t think it will help any due to the way G_plot_point()
works.

---

Also, there is still no one-step way to change a cell value based on
its xy coordinates.

So you want a non-interactive d.rast.edit, like v.edit?

A couple of ways to do that, 1) get input x1,y1,cat1 … from stdin,
qsort() based on y then x, read each line and if there’s an updated
cell in it replace it before writing the row back out. (as r.patch)
2) calculate the array column,row which contains each x,y and update
each in the order they arrive (as r.in.xyz). Probably the fastest way is
a combination of the two. Sort by y, and for each row update the column
values in the order they arrive, then write out the finished line and
move on to the next row.

However, if there is interest in doing agent and cellular modeling
within a GRASS platform, it will become necessary to be able to do
this to simulate dynamics. Patching may simply not be practical for
this kind of work. As Python becomes more used as a scripting
environment, this is an attractive kind of application to build in
Python, given its object oriented structure, and couple with GRASS.

The module probably isn’t so hard to write, in C. r.example should give
a nice start.

It would be interesting to rewrite r.example.c using the SWIG Python
interface, and see if that gains any simplicity over C, or if it just
adds complexity without simplicity – as you still need to go through
all the same chores.

Hamish

---

Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

---

grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

–
David Finlayson, Ph.D.
Operational Geologist

U.S. Geological Survey
Pacific Science Center
400 Natural Bridges Drive
Santa Cruz, CA 95060, USA

Tel: 831-427-4757, Fax: 831-427-4748, E-mail: dfinlayson@usgs.gov

[cmbarton]

Hi David,

Thanks for the ideas.

On 6/21/07 10:51 AM, "David Finlayson" <dfinlayson@usgs.gov> wrote:

Michael -

I often couple GRASS to numerical models or complex grid-based computations
using the following method:

Step 1: Export the environmental grids needed for the model using r.out.arc
(an easy format to parse in Python) or a simple binary format.
Step 2: Read the data into a Python class that represents the grid as a 2D
array (numeric is great for this). The class also has a few helper methods
such as: read_grid(), write_grid(), get(x, y), set(x,y), etc.
Step 3: Run the model
Step 4: Dump the python grid back to the ASCII Grid format (or raw binary)
Step 5: Read the data back into GRASS using r.in.arc (etc)

In my Mediterranean Landscape Dynamics project, we are trying a coupled
model approach. Our approach is have a OO-based ABM (in Java in this case)
do the decision modeling, and handle individual agents and their behavior,
but have C-based GRASS do the GIS/grid operations. As you point out below, C
is better at this--and a GIS like GRASS is optimized to handle sophisticated
processing of very large grids (in the 100K to millions of cells).

The idea is to let each platform do what it does best, rather than shoehorn
all the operations in to a single platform. Currently, it looks very
promising (I should see a new version this week or the next) and I'm
thinking that the same thing could be done with Python.

Michael

Most of the work in the model occurs in Steps 2-4 where there is no overhead
on the model due to carrying around the GIS baggage. You could use python, R,
Matlab, whatever.

As an unsolicited aside, C is actually a pretty good language for grid-based
models when one uses a numerical library such as the GNU Scientific Library to
handle the matrix (grid) stuff. Two bonuses to using C are the increase in
performance of the model and by the end of the project, you have learned
enough C to read and start using the GRASS source code directly. Not only
that, but there are lots of GIS-related libraries ( i.e., GDAL, Proj4,
Shapelib) that are easy to use from C once you learn the basics.

I resisted C for years because Python was good enough and C was "hard". Today
I think that was being penny-wise and pound foolish. There are many projects
in GIS where speed really matters. Whether you are filtering a very large
Lidar data set, or running a very small model millions of times (such as in a
Monte Carlo simulation), speed really matters in these cases. The huge library
of existing C code is real boon, too. I find that for small projects (< 10k
lines of code) C is a pretty good language; its warts don't show up until
project get really large. And like I said, it's nice to be able to customize
GRASS source when you need to.

Cheers,

David

On 6/14/07, Michael Barton <michael.barton@asu.edu> wrote:

Thanks very much. I'll need to digest all this useful information later
today, but it is helpful. I'll check r.in.xyz and your update to r.in.poly
and get back to you.

Michael

On 6/14/07 2:45 AM, "Hamish" <hamish_nospam@yahoo.com> wrote:

Michael Barton wrote:

We have finally had an initial production test of methods to create a
raster point using xy coordinates. If you remember, I've periodically
whined a bit about the lack of a way to create and change a raster
cell based on its xy coordinates instead of its cat value.

The 2 suggestions I got to work around this were to...

1) use r.in.poly to create the raster map; make an input file such
that for each point, the input is set to line (L) with 2 points that
are identical.

To make that method easier I have just added a "P" instruction to
r.in.poly in 6.3cvs. That uses G_plot_point(), which draws a 0 length
line using the same method, so probably it doesn't fix your 2 cell if
on bound problem (?? try it and see).

2) use v.in.ascii and then use v.to.rast to create the raster map.

Note the v.to.rast code is ~95% a clone of the r.in.poly code, so results
(and problems) should be mostly the same. Do you see the same multiple
cell trouble if the point falls on a bound using this method?

It probably needs to be noted in the v.to.rast help page if it does, at
minimum.

3) use r.in.xyz (I expect this will be the fastest method)

4) r.mapcalc outmap='if( (abs(x() - $x)) < ewres()) \
      && (abs(y() - $y) < nsres(), $val, null() )'
[or something like that, for a single point; yuk]

5) do the gridding yourself and create a r.in.ascii file.

all these would need to be followed by a call to r.patch.
(see d.rast.edit.tcl method)

We are coupling an agent-based simulation to GRASS, such that the
simulated land use activities will affect land cover, which in turn
will affect a landscape evolution routine.

are you updating a sequential string of cells, blocks of cells, or
random cells at each update?

Optimizing the speed of this is critical, and so we prefer to reduce
the number of steps needed to accomplish this. This means that the
r.in.poly approach is the one we prefer. However, in tests conducted
today, we ran into an unexpected problem.

If the coordinate of the point to be created happens to fall on a cell
boundary, we get a 2-cell 'line' instead of a point even though the
beginning and end points are identical. Apparently, r.in.poly looks
for a cell center closest to end point one and then searches for a
cell center for end point 2 that is different from the first cell,
resulting in a 2 cell line. This does not happen if the point does not
fall on a cell boundary. Attached is the raster file created from the
following r.in.poly input. The first input creates the 2 yellow cells,
the other 2 create the red and green single cells.

r.in.poly/raster.c defines the "cont" continue line fn used by
G_setup_plot() (and thus all following G_plot_*() commands) with
G_bresenham_line(). That in turn uses the "dot" fn to draw the points;
which "dot" is used depends on the result of the getformat() scan of
the cat values (& I'm about foggy about what that's all about).

for info on the Bresenham line rasterization method see:
  Bresenham's line algorithm - Wikipedia
<http://en.wikipedia.org/wiki/Bresenham&#39;s_line_algorithm&gt;

This is a potential problem, since an agent might not know that the
the coordinates it chooses for a field to farm or graze like on a cell
boundary. It is OK if r.in.poly substitutes the nearest grid cell
center to make the raster point for the xy coordinates that the agent
chooses, but not OK if it picks 2 different cells instead of one.

Try r.in.xyz, it shouldn't suffer from this: on one side of the cell it
tests for <= and on the other side <. So a point on a bound always fall
in the <= cell. The exception is the outer E,W raster boundaries which
will make values fall "into" the raster if on the bound.

r.in.xyz doesn't let you add cat labels, but you can do that manually by
creating a cats/ file. See spearfish60/PERMANENT/cats/landuse for an
example, it's very simple.

The simplest solution is to allow r.in.poly to make points as well as
lines and areas. I'm not sure why it doesn't allow points in the first
place. It seems like an oversight. Is this difficult to implement?

done, but I don't think it will help any due to the way G_plot_point()
works.

* * * *

Also, there is still no one-step way to change a cell value based on
its xy coordinates.

So you want a non-interactive d.rast.edit, like v.edit?

A couple of ways to do that, 1) get input x1,y1,cat1 ... from stdin,
qsort() based on y then x, read each line and if there's an updated
cell in it replace it before writing the row back out. (as r.patch)
2) calculate the array column,row which contains each x,y and update
each in the order they arrive (as r.in.xyz). Probably the fastest way is
a combination of the two. Sort by y, and for each row update the column
values in the order they arrive, then write out the finished line and
move on to the next row.

However, if there is interest in doing agent and cellular modeling
within a GRASS platform, it will become necessary to be able to do
this to simulate dynamics. Patching may simply not be practical for
this kind of work. As Python becomes more used as a scripting
environment, this is an attractive kind of application to build in
Python, given its object oriented structure, and couple with GRASS.

The module probably isn't so hard to write, in C. r.example should give
a nice start.

It would be interesting to rewrite r.example.c using the SWIG Python
interface, and see if that gains any simplicity over C, or if it just
adds complexity without simplicity -- as you still need to go through
all the same chores.

Hamish

__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

__________________________________________
Michael Barton, Professor of Anthropology and Director of Graduate Studies
School of Human Evolution & Social Change
Center for Social Dynamics and Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton