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

[Agustin Lobo]

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.

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


-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster