[R-sig-Geo] overlay polygons (from shp) and raster (from geotif)

[Roger Bivand]

On Tue, 23 Oct 2007, Agustin Lobo wrote:

> Hi there!
>
> I'm trying to do the overlay of polygons on raster.
>
> I do:
>
> pols <- readOGR("../AllTransectPolygons02", layer="AllTransectPolygons02")
> a <- readGDAL("t1s2comb/t1s2combfim95.tif") #geotif file
>
> str(a) shows the correct proj info:
>   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
>   .. .. ..@ projargs: chr " +proj=utm +zone=18 +south +ellps=WGS72
> +units=m +no_defs"
>
> My goal is to get, for each polygon in pols,
> the surface covered by each value of a.
>
> Then:
>
> delme <- overlay(pols, a)
> delme
>      AREA PERIMETER Site
> NA     NA        NA <NA>
> NA.1   NA        NA <NA>
>
> which is not too useful and
>
> delme2 <- overlay(a, pols)
>
> takes for ever and what I finally get is:
>
> > str(delme2)
>  num [1:207024] NA NA NA NA NA NA NA NA NA NA ...
> > delme2[!is.na(delme2)]
>  [1]  86  86  86  86  87  87  87  87  88  88  88  88  89  89  89  89
> 90  90
> [19]  90  90  90  91  91  92  92  93  93  93  93  94  94  94  95  95  95  96
> [37]  96  96  96  98 101 101 101 101 102 102 102 102 104 104 104 104 104 104
> [55] 104 104 105 105 105 105 105 105 106 106 106 106
>
> that is, just the ids of the polygons on the cells of a.

This is what you would need, using this column to make a 
SpatialPixelsDataFrame out of the input data, then doing a tapply or 
similar to tabulate the categories.

But with 200K cells, large cells, and small polygons, the point-in-polygon 
approach is converting the cells to their cell centre points, and missing 
most of the cell-polygon overlaps, that is those where the cell centre 
falls outside the polygon. Large cells and small polygons is not really 
the usual case. In addition, it might be helpful to subset pols first to 
avoid looping over polygons that cannot belong to a.

I'm not sure how starspan does polygon/cell overlay - if it measures the 
overlap area, you might find that more helpful. If not, reduce pols to the 
a area, and resample a to a finer resolution so that multiple cells fit 
inside each polygon (change g.region in GRASS, for example). It will take 
time to do, because it does a C point-in-polygon operation for each 
polygon (in nested lapply() calls).

Hope this helps,

Roger

>
> Note that the extent of pols is much larger than
> the one of a (i.e., may polygons are outside the raster and I'm not
> interested on those). Also, cells are much larger then the width of
> the polygons (I can send a jpg with an example).
>
> So it's clear I'm doing something wrong... any advice?
>
> Thanks
>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no