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