[R-sig-Geo] overlay polygons (from shp) and raster (from geotif)
Roger Bivand
Roger.Bivand at nhh.no
Wed Oct 24 10:17:44 CEST 2007
On Wed, 24 Oct 2007, Agustin Lobo wrote:
> Roger,
>
> A raster with the required resolution would be huge, although I can make an
> script setting the region for each polygon, doing the operation and
> combining all the tables afterwards. This is an alternative.
Look at something like this to get a list of bounding boxes:
lapply(slot(pols, "polygons"), function(x) sp:::bbox.Polygons(x))
which you could run out through g.region to set high resolution for each
Polygons object in turn.
>
> But I still think that a vector-based solution is more appropriate for this
> problem. I thought I could use joinPolys from PBSmapping, it has an
> intersection operator.
>
> My problem is that I've not found the way of making SpatialPolygons from
> SpatialPixels as you suggest.
as(as(a, "SpatialPixels"), "SpatialPolygons")
or equivalently
sp:::as.SpatialPolygons.SpatialPixels
but you need to watch the IDs.
> I've tried with adehabitat, but getcontour loses isolated pixels and
> does not transfer pixel values to the polygons. How could I convert the
> "SpatialGridDataFrame" into an SpatialPolygons object? Then I would need
> to convert to PBSmapping Polysets...
>
I would go straight to gpclib, which is what PBSmapping is using (in its
own form) internally.
For each Polygons object (slot(polys, "polygons")), find the raster cells
that intersect the bounding box of the object. Build gpclib polygon
rectangles by hand, and intersect with the Polygons object using code from
maptools, for example from checkPolygonsHoles(), to build the gpclib
objects. Then analyse the result to retrieve the intersections, possibly
by area.
> BTW, there is no vector intersection in grass either, am I wrong?
>
I think that there is v.overlay in 6.2 and 6.3 CVS, which might work from
r.to.vect (without corner smoothing), but I do not know how you will be
able to count the raster categories.
Roger
> Thanks,
>
> Agus
>
> Roger Bivand escribió:
>> 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
More information about the R-sig-Geo
mailing list