Skip to content

problem with rasterToPolygons x worldclim

6 messages · Karla Shikev, Michael Sumner

Karla Shikev
# 
Dear all,

Here's an issue that is related to the previous one. In the commands below
I'm trying to make a polygon for all the regions within a range of mean
annual temps. However, rasterToPolygons will draw a square around each
value in the original raster, rather than providing me the actual polygon,
given that resolution. Any hints? Again, any help will be greatly
appreciated.

Karla

_____

library(raster)

r<-getData('worldclim', var='bio', res=10)$bio1

e <- extent(-140,-100, 50, 60)

xx<-crop(r,e)

slice1<-rasterToPolygons(xx, fun = function(x) {x > 40 & x < 60})

plot(slice1)
Michael Sumner
# 
Read the help for rasterToPolygons, alternatively try chaining
rasterToContour and rgeos::gPolygonize.

The latter may need coaxing if your edges meet the raster extents.

Cheers, Mike
On Tue, 5 Apr 2016, 05:32 Karla Shikev <karlashikev at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      Karla Shikev
    

    

      
      #
    

  


  

      
Hi Michael,

Thanks for the tips. I tried the help for rasterToPolygons, but none of the
options (e.g. dissolve=TRUE) made any difference. I tried gPolygonize and
it worked, except for (as you predicted), if the edges meet the raster
extents - but that definitely was an advance!

I really appreciate your help and if you just point the way I can try to go
after what needs to be done, but right now I?m stumped.

best, Karla.

        
On Mon, Apr 4, 2016 at 5:38 PM, Michael Sumner <mdsumner at gmail.com> wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
On Wed, 6 Apr 2016 at 23:45 Karla Shikev <karlashikev at gmail.com> wrote:

        

            

      
      
      
    

        
I see the issue, the fun argument to rasterToPolygons is masking out values
from that interval, but it's not setting them all to the same value - so
you still get individual pixel polygons unless you mask on presence in the
interval or not (as a binary):

library(raster)

r<-getData('worldclim', var='bio', res=10)$bio1

e <- extent(-140,-100, 50, 60)

xx <- crop(r,e)
## mask out pixel values first
xx <- xx > 40 & xx < 60

## 576 polygons
slice1<- rasterToPolygons(xx, fun = function(x) {x == 1})
## 1 polygon
slice2 <- rasterToPolygons(xx,  fun = function(x) {x == 1}, dissolve = TRUE)


Still it's not a very nice polygon, if I can I'll try a different way.

Cheers, Mike.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

  
  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
On Thu, 7 Apr 2016 at 00:22 Michael Sumner <mdsumner at gmail.com> wrote:

        

            

      
      
      
    

        
You can use tricks to fill the outer edge so you get a valid contour line
all the way around (it's not always going to work - see here for the
explanation by whuber:
http://gis.stackexchange.com/questions/61550/colouring-areas-between-vector-contours
)

Is this better? Maybe - there are lots of other ways - maybe a buffer on
this gives roughly what you need. Maybe you only need a total area so
summing pixels after filtering is better?

Would be nice to have contouring that reliably produced polygons, but it's
a tricky problem. Other kinds of shape-finding might be better, with
alphahull on the pixel points perhaps, but you quickly get into "geometry
fudger" territory here rather than objective measures.

HTH

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

  
  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
On Thu, 7 Apr 2016 at 00:33 Michael Sumner <mdsumner at gmail.com> wrote:

        

            

      
      
      
    

        
Sorry, forgot the code!

library(raster)

r<-getData('worldclim', var='bio', res=10)$bio1

e <- extent(-140,-100, 50, 60)

xx <- crop(r,e)

## buffer out a little
xx2 <- extend(xx, extent(xx) + res(xx) * 4, value = NA_real_)

## set missing to less than interval
xx2[is.na(xx2)] <- cellStats(xx2, min)
## now contour
cl <- rasterToContour(xx2 > 40 & xx2 < 60, level = 1)

pp <- rgeos::gUnionCascaded(rgeos::gPolygonize(cl))
pp
plot(pp, col = "grey")

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia