sampling a raster within a polygon

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

sampling a raster within a polygon

Andy Bunn-3
Hi all, I’ve followed the thread on the updates to raster and would like to acknowledge the wonderful work done by Robert and others in making raster. It is an amazing piece of software. Thanks to all who help people like me do my job using the fruits of your labors!

Anyways, I’m working with some largish rasters and running into a step that is more computationally intensive than I would like. I’m trying to randomly sample a raster within the bounds of a polygon. E.g., in the example below I’d like to apply sampleRandom() on the raster r but constrain the sample to points within the polygons defined by xy.sp.


    require(raster)
   
    xy = cbind(
      x = c(13.4, 13.4, 13.6, 13.6, 13.4),
      y = c(48.9, 49, 49, 48.9, 48.9)
    )
    hole.xy <- cbind(
      x = c(13.5, 13.5, 13.45, 13.45, 13.5),
      y = c(48.98, 48.92, 48.92, 48.98, 48.98)
    )
   
    xy.sp <- SpatialPolygons(list(
      Polygons(list(Polygon(xy),
                    Polygon(hole.xy, hole = TRUE)), "1"),
      Polygons(list(Polygon(xy + 0.2),
                    Polygon(xy + 0.35),
                    Polygon(hole.xy + 0.2, hole = TRUE)), "2")
    ))
   
    r <- raster(nrow=100, ncol=100, ext=extent(xy.sp), resolution=0.01)
    r[] <- runif(ncell(r))
   
    plot(xy.sp,col="grey")
    plot(r,add=T,alpha=0.5)



I can accomplish this with a mask  via:
    sampleRandom(mask(r, xy.sp),size = 10)
   
However, I have many different polygons to apply over a big raster and the mask() is taking a long time. Is there a better way to go about this?

Thanks in advance, A    


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: sampling a raster within a polygon

obrlsoilau
Hi Andy,

I've gotten around this before using package 'sf' and list-columns,
here's a workflow using your data:

library(sf)
library(tidyverse)

xy.sf <- st_as_sf(xy.sp) %>%
  mutate(ID = 1:nrow(.)) # need at least one attribute

xy.sf.s <- xy.sf %>%
  # get intersecting cells
  mutate(intersecting_cells =
           # this is my go-to for row-wise functions atm, open to improvements
           split(., 1:nrow(.)) %>%
           map(., function(g) {
             poly_sp <- as(g, 'Spatial') # can we prioritise raster
support for sf pls :P
             cells   <- as.integer(unlist(cellFromPolygon(r, poly_sp)))
             cells   <- if (length(cells) == 0) { NA } else { cells }
             }) %>%
           setNames(., NULL)) %>%
  filter(!(is.na(intersecting_cells))) %>%
  # drop geometry here - its easier. output is just a df with several list cols
  st_set_geometry(., NULL) %>%
  # randomly subsample each polygon's cells
  mutate(sample =
           split(., 1:nrow(.)) %>%
           # you could store desired sample n in another column, but
lets go with 10:
           map(function(row) {
             sample(unlist(row[ , 'intersecting_cells']), size = 10,
replace = FALSE)
           })) %>%
  # get values of those cells
  mutate(values =
          split(., 1:nrow(.)) %>%
          map(function(row) {
            raster::extract(r, y = unlist(row[ , 'sample'])) %>%
              setNames(., NULL)
          }))

# proof that samples are where they should be
samp_1 <- xyFromCell(r, xy.sf.s$sample[[1]]) %>%
  st_multipoint(.) %>%
  st_sfc()
samp_2 <- xyFromCell(r, xy.sf.s$sample[[2]]) %>%
  st_multipoint(.) %>%
  st_sfc()
plot(samp_1, add = T, col = 'red', pch = 19)
plot(samp_2, add = T, col = 'blue', pch = 19)

HTH,
Lauren

On Tue, Oct 24, 2017 at 2:57 AM, Andy Bunn <[hidden email]> wrote:

> Hi all, I’ve followed the thread on the updates to raster and would like to acknowledge the wonderful work done by Robert and others in making raster. It is an amazing piece of software. Thanks to all who help people like me do my job using the fruits of your labors!
>
> Anyways, I’m working with some largish rasters and running into a step that is more computationally intensive than I would like. I’m trying to randomly sample a raster within the bounds of a polygon. E.g., in the example below I’d like to apply sampleRandom() on the raster r but constrain the sample to points within the polygons defined by xy.sp.
>
>
>     require(raster)
>
>     xy = cbind(
>       x = c(13.4, 13.4, 13.6, 13.6, 13.4),
>       y = c(48.9, 49, 49, 48.9, 48.9)
>     )
>     hole.xy <- cbind(
>       x = c(13.5, 13.5, 13.45, 13.45, 13.5),
>       y = c(48.98, 48.92, 48.92, 48.98, 48.98)
>     )
>
>     xy.sp <- SpatialPolygons(list(
>       Polygons(list(Polygon(xy),
>                     Polygon(hole.xy, hole = TRUE)), "1"),
>       Polygons(list(Polygon(xy + 0.2),
>                     Polygon(xy + 0.35),
>                     Polygon(hole.xy + 0.2, hole = TRUE)), "2")
>     ))
>
>     r <- raster(nrow=100, ncol=100, ext=extent(xy.sp), resolution=0.01)
>     r[] <- runif(ncell(r))
>
>     plot(xy.sp,col="grey")
>     plot(r,add=T,alpha=0.5)
>
>
>
> I can accomplish this with a mask  via:
>     sampleRandom(mask(r, xy.sp),size = 10)
>
> However, I have many different polygons to apply over a big raster and the mask() is taking a long time. Is there a better way to go about this?
>
> Thanks in advance, A
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: sampling a raster within a polygon

Bacou, Melanie
In reply to this post by Andy Bunn-3

Andy,
Simple use `extract()` instead of `mask()`, and then randomly sample cells in each polygon. Should be much faster, e.g.:

xy.r <- extract(r, xy.sp)
xy.r.sample <- lapply(xy.r, sample.int, n=10)

--Mel.

On 10/23/2017 12:57 PM, Andy Bunn wrote:
Hi all, I’ve followed the thread on the updates to raster and would like to acknowledge the wonderful work done by Robert and others in making raster. It is an amazing piece of software. Thanks to all who help people like me do my job using the fruits of your labors!

Anyways, I’m working with some largish rasters and running into a step that is more computationally intensive than I would like. I’m trying to randomly sample a raster within the bounds of a polygon. E.g., in the example below I’d like to apply sampleRandom() on the raster r but constrain the sample to points within the polygons defined by xy.sp.


    require(raster)
    
    xy = cbind(
      x = c(13.4, 13.4, 13.6, 13.6, 13.4),
      y = c(48.9, 49, 49, 48.9, 48.9)
    )
    hole.xy <- cbind(
      x = c(13.5, 13.5, 13.45, 13.45, 13.5),
      y = c(48.98, 48.92, 48.92, 48.98, 48.98)
    )
    
    xy.sp <- SpatialPolygons(list(
      Polygons(list(Polygon(xy),
                    Polygon(hole.xy, hole = TRUE)), "1"),
      Polygons(list(Polygon(xy + 0.2),
                    Polygon(xy + 0.35),
                    Polygon(hole.xy + 0.2, hole = TRUE)), "2")
    ))
    
    r <- raster(nrow=100, ncol=100, ext=extent(xy.sp), resolution=0.01)
    r[] <- runif(ncell(r))
    
    plot(xy.sp,col="grey")
    plot(r,add=T,alpha=0.5)



I can accomplish this with a mask  via: 
    sampleRandom(mask(r, xy.sp),size = 10)
    
However, I have many different polygons to apply over a big raster and the mask() is taking a long time. Is there a better way to go about this?

Thanks in advance, A     


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: sampling a raster within a polygon

Bacou, Melanie

Sorry the exact code is:

xy.r <- extract(r, xy.sp)
xy.r.sample <- lapply(xy.r, sample, 10)


On 10/24/2017 01:34 PM, Bacou, Melanie wrote:

Andy,
Simple use `extract()` instead of `mask()`, and then randomly sample cells in each polygon. Should be much faster, e.g.:

xy.r <- extract(r, xy.sp)
xy.r.sample <- lapply(xy.r, sample.int, n=10)

--Mel.

On 10/23/2017 12:57 PM, Andy Bunn wrote:
Hi all, I’ve followed the thread on the updates to raster and would like to acknowledge the wonderful work done by Robert and others in making raster. It is an amazing piece of software. Thanks to all who help people like me do my job using the fruits of your labors!

Anyways, I’m working with some largish rasters and running into a step that is more computationally intensive than I would like. I’m trying to randomly sample a raster within the bounds of a polygon. E.g., in the example below I’d like to apply sampleRandom() on the raster r but constrain the sample to points within the polygons defined by xy.sp.


    require(raster)
    
    xy = cbind(
      x = c(13.4, 13.4, 13.6, 13.6, 13.4),
      y = c(48.9, 49, 49, 48.9, 48.9)
    )
    hole.xy <- cbind(
      x = c(13.5, 13.5, 13.45, 13.45, 13.5),
      y = c(48.98, 48.92, 48.92, 48.98, 48.98)
    )
    
    xy.sp <- SpatialPolygons(list(
      Polygons(list(Polygon(xy),
                    Polygon(hole.xy, hole = TRUE)), "1"),
      Polygons(list(Polygon(xy + 0.2),
                    Polygon(xy + 0.35),
                    Polygon(hole.xy + 0.2, hole = TRUE)), "2")
    ))
    
    r <- raster(nrow=100, ncol=100, ext=extent(xy.sp), resolution=0.01)
    r[] <- runif(ncell(r))
    
    plot(xy.sp,col="grey")
    plot(r,add=T,alpha=0.5)



I can accomplish this with a mask  via: 
    sampleRandom(mask(r, xy.sp),size = 10)
    
However, I have many different polygons to apply over a big raster and the mask() is taking a long time. Is there a better way to go about this?

Thanks in advance, A     


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: sampling a raster within a polygon

Andy Bunn-3
Thanks Mel and Lauren for your replies. Very helpful indeed.

From: "Bacou, Melanie" <[hidden email]>
Reply-To: "[hidden email]" <[hidden email]>
Date: Tuesday, October 24, 2017 at 10:37 AM
To: Andy Bunn <[hidden email]>, R-sig-Geo <[hidden email]>
Subject: Re: [R-sig-Geo] sampling a raster within a polygon


Sorry the exact code is:

xy.r <- extract(r, xy.sp)
xy.r.sample <- lapply(xy.r, sample, 10)

On 10/24/2017 01:34 PM, Bacou, Melanie wrote:

Andy,
Simple use `extract()` instead of `mask()`, and then randomly sample cells in each polygon. Should be much faster, e.g.:

xy.r <- extract(r, xy.sp)
xy.r.sample <- lapply(xy.r, sample.int, n=10)
--Mel.
On 10/23/2017 12:57 PM, Andy Bunn wrote:

Hi all, I’ve followed the thread on the updates to raster and would like to acknowledge the wonderful work done by Robert and others in making raster. It is an amazing piece of software. Thanks to all who help people like me do my job using the fruits of your labors!



Anyways, I’m working with some largish rasters and running into a step that is more computationally intensive than I would like. I’m trying to randomly sample a raster within the bounds of a polygon. E.g., in the example below I’d like to apply sampleRandom() on the raster r but constrain the sample to points within the polygons defined by xy.sp.





    require(raster)



    xy = cbind(

      x = c(13.4, 13.4, 13.6, 13.6, 13.4),

      y = c(48.9, 49, 49, 48.9, 48.9)

    )

    hole.xy <- cbind(

      x = c(13.5, 13.5, 13.45, 13.45, 13.5),

      y = c(48.98, 48.92, 48.92, 48.98, 48.98)

    )



    xy.sp <- SpatialPolygons(list(

      Polygons(list(Polygon(xy),

                    Polygon(hole.xy, hole = TRUE)), "1"),

      Polygons(list(Polygon(xy + 0.2),

                    Polygon(xy + 0.35),

                    Polygon(hole.xy + 0.2, hole = TRUE)), "2")

    ))



    r <- raster(nrow=100, ncol=100, ext=extent(xy.sp), resolution=0.01)

    r[] <- runif(ncell(r))



    plot(xy.sp,col="grey")

    plot(r,add=T,alpha=0.5)







I can accomplish this with a mask  via:

    sampleRandom(mask(r, xy.sp),size = 10)



However, I have many different polygons to apply over a big raster and the mask() is taking a long time. Is there a better way to go about this?



Thanks in advance, A





_______________________________________________

R-sig-Geo mailing list

[hidden email]<mailto:[hidden email]>

https://stat.ethz.ch/mailman/listinfo/r-sig-geo




        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo