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

[Roger Bivand]

On Tue, 23 Oct 2007, Agustin Lobo wrote:

> Thanks,
>
> but note that:
>
>>  "Large cells and small
>>  polygons is not really the usual case."
>
> ...except when you carry out ground work for satellite imagery!
>
> According to the doc, starspan counts pixels intersected by the vector
> provided the intersection is larger than a user-define threshold, thus would 
> not solve the problem either (but have not actually tried it yet).
>
> In R, may be I could
> 1. Vectorize the raster cells (only those that are not NA in delme), thus 
> creating av.

There are facilities in sp for this - to make SpatialPolygons from 
SpatialPixels.

> 2. Intersect the 2 polygons, thus creating a new set of polygons (pols2)

This would possibly involve using gpclib directly - there is code in 
maptools for going between SpatialPolygons and the gpclib representation, 
but it is still inside functions.

> 3. Use overlay with the new polygon pols2 and the vectorized version (av) of 
> the raster a.

There is no polygon/polygon overlay method in sp - you'd need to keep 
track of which polygon is which in step 2.

I think I still prefer resampling a in a GIS to a higher resolution and 
using what exists.

Roger

>
> What do you think?
>
> Agus
>
>
> Roger Bivand escribió:
>>  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