0

I want to crop a raster and then save each peace into a list (or another way) with the pieces of the raster.

The idea:

  • To preserve the covered area as equal for all of the pieces.
  • To be possible for the user to set a percentage of extent area to be cropped. I.e. 10, 20, 30% of the all-area should be the extent of the grid.

A piece of reproducible code can be seen below: It has been developed by Chris and Robert Hijmans at this question here.

library(terra)
r <- rast(volcano)
n <- 10

Get the starting cells of interest

rows <- seq(1, nrow(r), by=n)
cols <- seq(1, ncol(r), by=n)    
cells <- cellFromRowColCombine(r, rows, cols)

Get the coordinates

upper-left coordinates of the starting cells

xy <- xyFromCell(r, cells)
rs <- res(r)
xy[,1] <- xy[,1] - rs[1]/2
xy[,2] <- xy[,2] + rs[2]/2

add the lower-right coordinates of the end cell

xy <- cbind(xy[,1], xy[,1] + n*rs[1], xy[,2] - n*rs[2], xy[,2])

And loop

x <- lapply(1:nrow(xy), function(i) {
         crop(r, xy[i,])
       })

Verify

e <- lapply(x, \(i) ext(i) |> as.polygons()) |> vect()
plot(r)
lines(e, col="blue", lwd=2)

enter image description here

Note the right side of the cropped raster. It seems bad to perform statistics on these polygons.

Also: I would like to:

  • the n on the beginning meaning a percentage of the area to be cropped in little rasters
  • The user would enter 10, 20, or 30 as a number, for instance, meaning a percentage on the command to decide the size of the "window" to be cropped
  • To preserve an equal space on the window, and to do that the command should vary a little between 10 and 20, for instance, to preserve the symmetry.

1 Answers1

0

You should really try to focus on the question at hand, and show what you have tried. Because it remains hard to know what you want exactly. Anyway, with n=20% and a matrix of 40 (columns) by 60 (rows)

cols = 40
rows = 60
n = 20 

You can find block size like this:

h = sqrt(n/100) * rows
w = sqrt(n/100) * cols 

"proof"

p = h * w
100 * p / (cols * rows) 
# 20

Of course h and w are not integers

> h
#[1] 26.83282
> w
#[1] 17.88854

And you would get

c(seq(1, cols, w) |> round(), cols)
#[1]  1 19 37 40
c(seq(1, rows, h) |> round(), rows)
[1]  1 28 55 60

Another approach. First compute the number of blocks. With 20% that is five. In most cases, that should really be an even number, unless you take entire rows, or columns, so you cannot get exactly that.

# number of blocks
nb <- 100/n 
# block rows and cols
brow <- round(nb * rows/(rows+cols))
brow
#[1] 3
bcol <- round(nb / brow)
bcol
#[1] 2

And now you can perhaps use:

round(c(seq(1, rows, rows/brow), rows))
#[1]  1 21 41 60
round(c(seq(1, cols, cols/bcol), cols))
#[1]  1 21 40
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • `S Pavoine(2008) Gower's distance allows an efficient treatment of missing data and the inclusion of variable...` suggests setting 'excess' areas of a 'user's ROI' to NA once the calc'd impact of user's c(0.1,0.2,0.3) is evaluated and let Gower's treatment take over. The language suggest a trajectory toward a service portal whose viability is mediated by the data, user options, how those options are evaluated. – Chris Apr 09 '22 at 23:38