[GRASS-user] Calculating a raster "polygon density" map

[Matthew Perry]

Brent,

On 11/20/06, Brent Wood <b.wood at 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 me :-) using a XYZ
> dataset of grid centroids & whatever Z value I want to use
>

Thanks 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