Can I suggest that you post this kind of 'tip' on the grass wiki, either in FAQs or in the section on how to do things (I forget the title). It's much easier to track them down later.
Thanks, Jerry
---- Original message ----
Date: Thu, 4 Jan 2007 22:43:44 -0800
From: "Matthew Perry" <perrygeo@gmail.com>
Subject: Re: [GRASS-user] Calculating a raster "polygon density" map
To: b.wood@niwa.co.nz
Cc: grassuser@grass.itc.itBrent,
On 11/20/06, Brent Wood <b.wood@niwa.co.nz> wrote:
Matthew Perry wrote:
> Hey folks,
>
> I'm trying to rasterize a polygon vector layer as a polygon density
> map. In other words instead of each cell containing NULL or the
> underlying vector attribute (v.to.rast), I'd like each raster cell to
> represent the percent of the cell's area covered by polygons.I've recently done a similar exercise with polygons representing the
swept area of commercial fishing trawls and a grid polygon layer.The "raster" layer is a table of square polygons forming a grid covering
the area of interest.A simple SQL joins the two tables (400,000 cells & 1,000,000 tow
polygons), giving 5,000,000 records & the intersections of the two
layers, ie: each tow is divided up into the cells it crosses. Given the
output includes cell_id, SQL queries to determine the aggregate swept
area within a cell, which can be many times the cell area, the area of
each cell which has or has not been fished, etc are trivial & pretty fast.So my solution was not a raster one as you suggest , but used a vector
layer of cells which functioned in a similar way. It is then trivial to
get the data into GRASS (for you) or GMT & R (for meusing a XYZ
dataset of grid centroids & whatever Z value I want to useThanks for the tip .. I knew I couldn't be the only one doing this!
I gather you're using postgis to do the sql joins. That idea crossed
my mind but I was hesitant to try to convert a huge raster into a
layer with many millions of polygons. (our raster has 933 million
cells)I attempted to do something similar using ogr by scanning each raster
cell, making a temporary polygon, doing the intersection to calculate
coverage area, then write the percent areal coverage to the output
raster cell. It seems to work for smaller rasters (slowly) but it
chocked on the huge raster datasets so we put the project on permanent
haitus.Anyways, I thought someone might find it useful as a (rough) starting
point for a "polygon density" utility :
http://perrygeo.googlecode.com/svn/trunk/gis-bin/poly_density.py--
Matthew T. Perry
http://www.perrygeo.net_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser
Gerald Nelson
Professor, Dept. of Agricultural and Consumer Economics
University of Illinois, Urbana-Champaign
office: 217-333-6465
cell: 217-390-7888
315 Mumford Hall
1301 W. Gregory
Urbana, IL 61801