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

[Michael Barton]

Hi David,

Thanks for the ideas.


On 6/21/07 10:51 AM, "David Finlayson" <dfinlayson at 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 at 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 at 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
>>> <http://en.wikipedia.org/wiki/Bresenham&#39;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 at 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