[GRASS-dev] Re: report on raster point creation

[Hamish]

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