Skip to content

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

6 messages · Roger Bivand, Agustin Lobo

Original
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
Roger Bivand
#
On Tue, 23 Oct 2007, Agustin Lobo wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
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?:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Tue, 23 Oct 2007, Agustin Lobo wrote:

        

            

      
      
      
    

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

            

      
      
      
    

        
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.

            

      
      
      
    

        
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

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
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.

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. 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...

BTW, there is no vector intersection in grass either, am I wrong?

Thanks,

Agus

Roger Bivand escribi?:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 24 Oct 2007, Agustin Lobo wrote:

        

            

      
      
      
    

        
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.

            

      
      
      
    

        
as(as(a, "SpatialPixels"), "SpatialPolygons")

or equivalently

sp:::as.SpatialPolygons.SpatialPixels

but you need to watch the IDs.

            

      
      
      
    

        
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.

            

      
      
      
    

        
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