[GRASS-user] strategy required

[Moritz Lennert]

On 22/06/08 13:23, John Field wrote:
> G'day Grass listers,
> 
> I'm very new at using Grass, and I'd be most appreciative if someone 
> could suggest a strategy to accomplish the following:
> 
> I have a map of part of Australia, subdivided into local government 
> areas (LGAs).  Both the coastline and the LGA boundaries are vector 
> layers, referenced by latitude and longitude.
> 
> I want to sample the map at each point on a grid (say at each 0.1 of lat 
> and 0.1 of long).  At a sample point I will place a rectangle of a given 
> size oriented in a NW-SE direction, with the westernmost corner of the 
> rectangle on the sample point.  I then want to calculate the proportion 
> of each LGA which is enclosed in the rectangle (there may be up to 6 or 
> so).
> 
> I'm not familiar enough with Grass to know how best to approach this.  
> I'd be very grateful for any suggested strategies.

One idea:

Preliminaries:

- check here: 
http://spreadsheets.google.com/pub?key=pB-ffIn4C2dTcTacV_W1cOA for the 
GUI paths to most commands mentioned below

- I'm not sure how reliable area calculations are in lat-long locations, 
maybe you could consider using a projected system
- you will need a SQL-compliant database backend for this, the easiest 
would probably be SQLITE as this is installed with GRASS, to do this use 
db.connect to change your default database connection and then g.copy 
your LGA layer, thus automatically creating an SQLite table 
(alternative: use db.copy to copy attribute table from default dbf to 
SQLite database and v.db.connect to reconnect your LGA layer to that new 
table)

And here's the "strategy":

- use v.db.addcol to add a column 'area' to your LGA layer
- use v.to.db to add the area of each LGA to the table
- create grid with v.mkgrid  - if you use the position=coor option, you 
can give an angle of rotation to your grid
- use v.overlay, operator=or to overlay grid and LGA; the result will be 
a layer with new entities created by the combination of all boundaries 
of both layers, and a table containing all fields of both attribute 
tables, prefixed a_ and b_ for each original layer respectively.
- use v.db.addcol (let's call that column 'newarea') and v.to.db (as 
above) to add the areas of the new entities to the attribute table of 
the new layer (make sure you use the same units as in the previous call 
to v.to.db)
- then use an SQL query to calculate the areas by grid cell (assuming 
that key columns are named 'cat' and that the grid is layer a and the 
LGA layer b, something like:

SELECT a_cat, b_cat, sum(newarea)/b_area as proportion_LGAarea from 
NameOfNewLayer group by a_cat, b_cat

(you need something like the sum(newarea) if the same LGA creates 
several new entities in one grid cell)

Moritz