overlay polygons (from shp) and raster (from geotif)

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.

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

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

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