3

The goal is to divide a raster file into quadrants in order to make some computation in each. I'd like to split the raster and create a raster brick with those parts. I don't need to get perfect clips of the raster so the cell with data in the border of the clip could be also present in the contiguous raster brick.

I'm thinking about creating a radial grid and then extract the raster by masking it with each of the polygons within. I found these interesting questions about the topic.

How can be done this in R ?

Here you have the layer (dem) to tryout.

enter image description here

  • Note 1: number of quadrants must be a parameter.
  • Note 2: raster brick must have same extension as the original raster.
  • Note 3: process must be done in R (a function would be perfect)

Update 10/06/2020: I hope this schema clarifies the purpose: A Function that given a raster and a number of sectors it creates a raster brick (with number of rasters = number of sectors) where each raster keeps just the values within the correspondent sector.

enter image description here

  • 1
    The polygons that define each quadrant only have three points (eg (0,0)->(xmax,0)->(xmax,ymax)) so you can define these polygons in eight very short lines. Have you tried? – Spacedman Jun 8 at 19:19
  • 1
    note: raster always have a rectangular extent. Therefore splitting by sectors is not so effiient. – radouxju Jun 14 at 19:58
  • 1
    Had you seen this: gis.stackexchange.com/questions/255025/… ? It is not a duplicate, but examples there must be useful. For each sector polygon i would use x:c(sin(k*2*pi/8)*width,0,sin((k+1)*2*pi/8)*width) y:c(cos(k*2*pi/8)*width,0,cos((k+1)*2*pi/8)*width), k ranges from 0 to 7, width is raster width. Those polygons will be surely bigger than needed, but must work just fine. – Javier JC Jun 15 at 19:33
  • 1
    You need to be much more clear on this. Is your input a square (NxN cell) raster stack and you want the output to be 8 smaller (N/2 x N/2) raster stacks which are square crops of the source and NA-masked by the circle segment? – Spacedman Jun 16 at 8:01
  • 1
    Or do you have a general raster stack, and given a centre point and radius you want the 8 cropped raster stacks with NA-masking outside the circle segment? – Spacedman Jun 16 at 8:03
1
+50

EDIT: Changed the function so it can now generate number of slices other than multiples of 4.

I've come up with a solution using the sf package and rotating points around the raster. Here's what a slice looks like.

Single Slice

It works perfectly, you can test it out yourself:

#' Sector Masking
#' @description Slices Raster object into sections
#'
#' @param raster Raster* object
#' @param segments Number of segments to be sliced
#'
#' @return RasterBrick containing slices in bands
#'
#' @import raster
#' @import sf
#' @import magrittr
#' 
sector_masking <- function(raster, segments){
  library(raster)
  library(sf)
  library(magrittr)
  
  # Get properties of raster
  dem <- brick(raster)
  ext <- extent(dem)
  crs <- as.character(crs(dem))
  
  # Calculate mid and max point of the raster
  half_x <- (ext@xmax - ext@xmin) / 2
  half_y <- (ext@ymax - ext@ymin) / 2
  middle <- st_point(c(ext@xmin+half_x, ext@ymin+half_y))
  max_ext <- st_point(c(ext@xmax, ext@ymax))
  
  # Rotates a starting point around the middle of the raster
  # These points are used to generate triangles which are masking the raster
  slices <- list()
  starting_point <- max_ext
  rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
  rotation <- segments/2
  for(i in 1:segments){
    next_point <- (starting_point - middle) * rot(pi/rotation) + middle
    slice_matrix <- rbind(middle,
                          starting_point,
                          next_point,
                          middle)
    slice_polygon <- slice_matrix %>% 
      list() %>% 
      st_polygon() %>% 
      st_sfc(crs = crs) %>% 
      st_sf
    slice_raster <- rasterize(slice_polygon, dem, mask = T)
    slices[[i]] <- slice_raster
    starting_point <- next_point
  }
  return(stack(slices))
}
| improve this answer | |
  • Thanks for the answer. I will try it as soon as possible. – César Arquero Jun 17 at 12:09
  • 1
    Does it work for you? – JonasV Jun 19 at 6:09
  • It works and the rasters keep extent, so that's perfect. I found that with less than 4 segments, the tool behaves weird, but I can try to figure it out. Good job. – César Arquero Jun 19 at 6:48
  • Oh yes, I didn't think too much about the edge cases. For three and two segments you would probably have to build the polygon in two steps. Otherwise parts of the raster will be cut off. – JonasV Jun 19 at 7:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.