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

Thu Jun 21 13:51:07 EDT 2007

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 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
> >
> >
> >> 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
>

-- 
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 at usgs.gov
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-dev/attachments/20070621/1f4da7f6/attachment.html