On Tue, 23 Oct 2007, Agustin Lobo wrote:
"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:
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:
num [1:207024] NA NA NA NA NA NA NA NA NA NA ...
[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
[37] 96 96 96 98 101 101 101 101 102 102 102 102 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