Michael -<br><br>I often couple GRASS to numerical models or complex grid-based computations using the following method:<br><br>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.
<br>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.<br>Step 3: Run the model
<br>Step 4: Dump the python grid back to the ASCII Grid format (or raw binary)<br>Step 5: Read the data back into GRASS using r.in.arc (etc)<br><br>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.
<br><br>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.<br><br>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.
<br><br>Cheers,<br><br>David<br><br><br><br><div><span class="gmail_quote">On 6/14/07, <b class="gmail_sendername">Michael Barton</b> <<a href="mailto:michael.barton@asu.edu">michael.barton@asu.edu</a>> wrote:</span>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Thanks very much. I'll need to digest all this useful information later<br>today, but it is helpful. I'll check
r.in.xyz and your update to r.in.poly<br>and get back to you.<br><br>Michael<br><br><br>On 6/14/07 2:45 AM, "Hamish" <<a href="mailto:hamish_nospam@yahoo.com">hamish_nospam@yahoo.com</a>> wrote:<br><br>> Michael Barton wrote:
<br>>> We have finally had an initial production test of methods to create a<br>>> raster point using xy coordinates. If you remember, I've periodically<br>>> whined a bit about the lack of a way to create and change a raster
<br>>> cell based on its xy coordinates instead of its cat value.<br>>><br>>> The 2 suggestions I got to work around this were to...<br>>><br>>> 1) use r.in.poly to create the raster map; make an input file such
<br>>> that for each point, the input is set to line (L) with 2 points that<br>>> are identical.<br>><br>> To make that method easier I have just added a "P" instruction to<br>> r.in.poly in
6.3cvs. That uses G_plot_point(), which draws a 0 length<br>> line using the same method, so probably it doesn't fix your 2 cell if<br>> on bound problem (?? try it and see).<br>><br>>> 2) use v.in.ascii
and then use v.to.rast to create the raster map.<br>><br>> Note the v.to.rast code is ~95% a clone of the r.in.poly code, so results<br>> (and problems) should be mostly the same. Do you see the same multiple<br>
> cell trouble if the point falls on a bound using this method?<br>><br>> It probably needs to be noted in the v.to.rast help page if it does, at<br>> minimum.<br>><br>><br>> 3) use r.in.xyz (I expect this will be the fastest method)
<br>><br>> 4) r.mapcalc outmap='if( (abs(x() - $x)) < ewres()) \<br>> && (abs(y() - $y) < nsres(), $val, null() )'<br>> [or something like that, for a single point; yuk]<br>><br>
> 5) do the gridding yourself and create a r.in.ascii file.<br>><br>> all these would need to be followed by a call to r.patch.<br>> (see d.rast.edit.tcl method)<br>><br>><br>>> We are coupling an agent-based simulation to GRASS, such that the
<br>>> simulated land use activities will affect land cover, which in turn<br>>> will affect a landscape evolution routine.<br>><br>> are you updating a sequential string of cells, blocks of cells, or<br>
> random cells at each update?<br>><br>><br>>> Optimizing the speed of this is critical, and so we prefer to reduce<br>>> the number of steps needed to accomplish this. This means that the<br>>>
r.in.poly approach is the one we prefer. However, in tests conducted<br>>> today, we ran into an unexpected problem.<br>>><br>>> If the coordinate of the point to be created happens to fall on a cell<br>
>> boundary, we get a 2-cell 'line' instead of a point even though the<br>>> beginning and end points are identical. Apparently, r.in.poly looks<br>>> for a cell center closest to end point one and then searches for a
<br>>> cell center for end point 2 that is different from the first cell,<br>>> resulting in a 2 cell line. This does not happen if the point does not<br>>> fall on a cell boundary. Attached is the raster file created from the
<br>>> following r.in.poly input. The first input creates the 2 yellow cells,<br>>> the other 2 create the red and green single cells.<br>><br>> r.in.poly/raster.c defines the "cont" continue line fn used by
<br>> G_setup_plot() (and thus all following G_plot_*() commands) with<br>> G_bresenham_line(). That in turn uses the "dot" fn to draw the points;<br>> which "dot" is used depends on the result of the getformat() scan of
<br>> the cat values (& I'm about foggy about what that's all about).<br>><br>> for info on the Bresenham line rasterization method see:<br>> <a href="http://en.wikipedia.org/wiki/Bresenham's_line_algorithm">
http://en.wikipedia.org/wiki/Bresenham's_line_algorithm</a><br>><br>><br>>> This is a potential problem, since an agent might not know that the<br>>> the coordinates it chooses for a field to farm or graze like on a cell
<br>>> boundary. It is OK if r.in.poly substitutes the nearest grid cell<br>>> center to make the raster point for the xy coordinates that the agent<br>>> chooses, but not OK if it picks 2 different cells instead of one.
<br>><br>> Try r.in.xyz, it shouldn't suffer from this: on one side of the cell it<br>> tests for <= and on the other side <. So a point on a bound always fall<br>> in the <= cell. The exception is the outer E,W raster boundaries which
<br>> will make values fall "into" the raster if on the bound.<br>><br>> r.in.xyz doesn't let you add cat labels, but you can do that manually by<br>> creating a cats/ file. See spearfish60/PERMANENT/cats/landuse for an
<br>> example, it's very simple.<br>><br>><br>>> The simplest solution is to allow r.in.poly to make points as well as<br>>> lines and areas. I'm not sure why it doesn't allow points in the first
<br>>> place. It seems like an oversight. Is this difficult to implement?<br>><br>> done, but I don't think it will help any due to the way G_plot_point()<br>> works.<br>><br>><br>>> * * * *
<br>>><br>>> Also, there is still no one-step way to change a cell value based on<br>>> its xy coordinates.<br>><br>> So you want a non-interactive d.rast.edit, like v.edit?<br>><br>> A couple of ways to do that, 1) get input x1,y1,cat1 ... from stdin,
<br>> qsort() based on y then x, read each line and if there's an updated<br>> cell in it replace it before writing the row back out. (as r.patch)<br>> 2) calculate the array column,row which contains each x,y and update
<br>> each in the order they arrive (as r.in.xyz). Probably the fastest way is<br>> a combination of the two. Sort by y, and for each row update the column<br>> values in the order they arrive, then write out the finished line and
<br>> move on to the next row.<br>><br>><br>>> However, if there is interest in doing agent and cellular modeling<br>>> within a GRASS platform, it will become necessary to be able to do<br>>> this to simulate dynamics. Patching may simply not be practical for
<br>>> this kind of work. As Python becomes more used as a scripting<br>>> environment, this is an attractive kind of application to build in<br>>> Python, given its object oriented structure, and couple with GRASS.
<br>><br>> The module probably isn't so hard to write, in C. r.example should give<br>> a nice start.<br>><br>> It would be interesting to rewrite r.example.c using the SWIG Python<br>> interface, and see if that gains any simplicity over C, or if it just
<br>> adds complexity without simplicity -- as you still need to go through<br>> all the same chores.<br>><br>><br>> Hamish<br><br>__________________________________________<br>Michael Barton, Professor of Anthropology
<br>School of Human Evolution & Social Change<br>Center for Social Dynamics & Complexity<br>Arizona State University<br><br>phone: 480-965-6213<br>fax: 480-965-7671<br>www: <a href="http://www.public.asu.edu/~cmbarton">
http://www.public.asu.edu/~cmbarton</a><br><br><br>_______________________________________________<br>grass-dev mailing list<br><a href="mailto:grass-dev@grass.itc.it">grass-dev@grass.itc.it</a><br><a href="http://grass.itc.it/mailman/listinfo/grass-dev">
http://grass.itc.it/mailman/listinfo/grass-dev</a><br></blockquote></div><br><br clear="all"><br>-- <br>David Finlayson, Ph.D.<br>Operational Geologist<br><br>U.S. Geological Survey<br>Pacific Science Center<br>400 Natural Bridges Drive
<br>Santa Cruz, CA 95060, USA<br><br>Tel: 831-427-4757, Fax: 831-427-4748, E-mail: <a href="mailto:dfinlayson@usgs.gov">dfinlayson@usgs.gov</a>