Skip to content

Polygons VS MultiLineString

7 messages · Antonio Rodriges, Michael Sumner, Frede Aakmann Tøgersen +1 more

Original
Antonio Rodriges
# 
Hello,

I use "rasterToContour" to derive isolines for raster. After that I
save them as GeoJSON (or ESRI Shapefile). However, it saves contours
as "MultiLineString"

Is it possible to save "contours" as "Polygons"?

Code:

to_c - is a raster
contours <- rasterToContour(to_c, levels = c(18, 33))
writeOGR(contours, "file.js", layer="", driver="GeoJSON")
Michael Sumner
# 
Try gPolygonize from rgeos, but it *really* depends on whether your lines
close cleanly.
On Wed, 4 Mar 2015 18:38 Antonio Rodriges <antonio.rrz at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Frede Aakmann Tøgersen
    

    

      
      #
    

  


  

      
Well, since the workhorse behind the rasterToContour() function is contourLines() from the lattice package don't expect the contour lines for a given level to form a single well defined polygon (the definition of which of course needs to be specified).

Try this:

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
x <- rasterToContour(r)
class(x)
str(x)

## the levels
levels(x at data$level)

plot(r)

plot(subset(x, level == 400), add=TRUE)

plot(subset(x, level == 200), add=TRUE)

Yours sincerely / Med venlig hilsen


Frede Aakmann T?gersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
<snip>

            

      
      
      
    

        
If your data don't build clean lines I would be inclined to classify the
raster data and then turn that into polygons.
But, that means the boundaries are not as smooth.

Using example data:

library(raster)
r <- raster(volcano)
cl <- rasterToContour(r, levels = c(110, 130))

## see how cl$level == 130 will polygonize, but not 110
library(rgeos)
p1 <- gPolygonize(subset(cl, level == 130))
p2 <- gPolygonize(subset(cl, level == 110))
plot(p1)  ## ok
is.null(p2) ## TRUE

And that is even before you decide if you want non-intersecting polygon
donuts etc.

So, forgetting contouring

## cut into our breaks
vals <- c(110, 130)
rc <- cut(r, c(vals, cellStats(r, max)))
plot(rc)
plot(cl, add = TRUE)

## another way to polygonize
rp <- rasterToPolygons(rc, dissolve = TRUE)
plot(rp, col = c("grey", "lightblue"))
rp$layer <- vals

## check all is well
plot(r)
plot(cl, add = TRUE)
contour(r, lev = vals, add = TRUE)
text(rp[2,], lab = rp$layer[2])

plot(rp, add = TRUE, lty = 2)

So there we end up with pretty rough boundaries, since the polygonizing
from cells is following pixel boundaries exactly. Maybe we can
smooth/disaggregate the raster, but it gets pretty expensive.

Another trick might be to buffer the raster with values outside the range,
and do contour individually
for each level.

So,

i <- 1
r1 <- r
r1[is.na(r1) | r1 < vals[i]] <- vals[i] - 10
r2 <- extend(r1, extent(r) + res(r) * 2, value = vals[i] - 10)

plot(r2)
plot(gPolygonize(cl1), add = TRUE, col = grey(0.5, alpha = 0.6))


Is that better, maybe - hard to generalize I imagine, and still you need to
clip your larger boundary with internal ones, and merge them into one data
set.

Thanks for the question though, it's given me some ideas to fix some things
I need to do. (Manifold, and I imagine other, GIS does this seamlessly so
it would be nice to have contour and contourLines upgraded to be able to
trace around the edges.)

Cheers, Mike.

        
On Wed, 4 Mar 2015 at 19:26 Frede Aakmann T?gersen <frtog at vestas.com> wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Antonio Rodriges
    

    

      
      #
    

  


  

      
Mike,

Thank you, that was a great example; I reproduced it.

However, if I am satisfied with the rasterToContour output, I think
there should be a way to directly convert it polygons. I am seeking
this method.

  
  


  


      

    

    
4 days later

    

      
      

        


  

        

      Antonio Rodriges
    

    

      
      #
    

  


  

      
The solution I've devised so far

contours <- rasterToContour(to_c, levels = c(18, 33))

res <- lapply(slot(contours, "lines"), function(x) lapply(slot(x,
"Lines"), function(y) slot(y, "coords")))
polies <- lapply(res, function (x) lapply(x, function (y) Polygon(y)))
pp_18 <- Polygons(polies[[1]], ID = "18")
pp_33 <- Polygons(polies[[2]], ID = "33")
all_pp<- SpatialPolygons(c(pp_18, pp_33))
df <- data.frame(value = c(18, 33), row.names = c("18", "33"))
SPDF <- SpatialPolygonsDataFrame(all_pp, df)

writeOGR(SPDF, jsfile, layer="", driver="GeoJSON")

This is in case you have only two values 18 and 33. If you have more,
create list and create Polygons also using lapply or for loop

2015-03-05 11:49 GMT+03:00 Antonio Rodriges <antonio.rrz at gmail.com>:

            

      
      
      
    

  
  


  


      

    

    
8 days later

    

      
      

        


  

        

      Mathieu Rajerison
    

    

      
      #
    

  


  

      
Hello,

My advice is : why don't you reclass your raster in categories according to
raster values intervals, then you polygonize your reclassified raster ?

It would be simpler and you would have a topologically correct layer.

That's the way I do it with GRASS but you all the functions also are
included in the raster package.

Mat



2015-03-09 12:44 GMT+01:00 Antonio Rodriges <antonio.rrz at gmail.com>: