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

César Arquero Cabral
  • 3,439
  • 1
  • 16
  • 53
  • 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 08 '20 at 19:19
  • Now that you said so it seems obvious...but I didn't. :-) I'll try that way. – César Arquero Cabral Jun 08 '20 at 21:30
  • 1
    note: raster always have a rectangular extent. Therefore splitting by sectors is not so effiient. – radouxju Jun 14 '20 at 19:58
  • Ler' say in other way: I want to put na values in all the sectors but the selected one. The raster brick should be like a cake splited into parts. – César Arquero Cabral Jun 14 '20 at 20:42
  • 1
    Had you seen this: https://gis.stackexchange.com/questions/255025/r-raster-masking-a-raster-by-polygon-also-remove-cells-partially-covered ? 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 '20 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 '20 at 08: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 '20 at 08:03
  • 1
    I've just did a simple schema of the goal. – César Arquero Cabral Jun 16 '20 at 08:38
  • 1
    What does it do if the number of "quadrants" is 7 or 5 or 6? – Spacedman Jun 16 '20 at 11:32
  • 1
    As a tease, I've written a function for the 4-sector case. Generalise to 8 sectors is easy. Not sure about any other value though. – Spacedman Jun 16 '20 at 11:49
  • 1
    Actually its quite easy to generalise to N sectors but you should specify where you want the sectors to start from. For example, starting from North and dividing the 360 degrees round into N equal sectors, but in your 4-sector case it starts from a direction that depends on the aspect ratio of the raster... – Spacedman Jun 16 '20 at 13:27
  • I think that that is the best way to do it: starting from north (0) and dividing by degrees clockwise. For example: a 7 sectors split should produce polygons at (0 --> 51,43 // 51,43 --> 102,86...) – César Arquero Cabral Jun 16 '20 at 13:54

1 Answers1

1

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))
}
JonasV
  • 3,754
  • 8
  • 24