Search code examples
rresolutiondimensionsrasterize

Issue with extracting rows from a Raster* in R

I'm running into an issue working with Raster* objects. I want to extract coordinates from the top and bottom rows of a raster, and use them to then create a new raster with values in only those cells. This is probably a bad explanation, so I'll provide an example:

r <- raster(xmn=-150, xmx=-120, ymx=60, ymn=30, ncol=36, nrow=18)
values(r) <- 1:ncell(r)

ncol_r <- ncol(r)
ncell_r <- ncell(r) 

N_nodes_cells_r <- 1:ncol_r                        # first row (N, for North)
S_nodes_cells_r <- ((ncell_r - ncol_r)+1):ncell_r  # final row (S, for South)

NS_nodes_r <- rbind(xyFromCell(r, cell=c(N_nodes_cells_r), spatial=T),
                    xyFromCell(r, cell=c(S_nodes_cells_r), spatial=T))

r_NS_nodes <- rasterize(x = NS_nodes_r, y = r)

plot(r_NS_nodes)

This works perfectly. But, I want to do it on a much larger scale. So, I replace the extent and ncol and nrow parameters with values equivalent to the real-world raster I'm working with. See below:

r <- raster(xmn=-117.2667, xmx=-97.13757, ymn=25.42458, ymx=33.16274, ncol=785, nrow=2042)
values(r) <- 1:ncell(r)

ncol_r <- ncol(r)
ncell_r <- ncell(r) 

N_nodes_cells_r <- 1:ncol_r
S_nodes_cells_r <- ((ncell_r-ncol_r)+1):ncell_r

NS_nodes_r <- rbind(xyFromCell(r, cell=c(N_nodes_cells_r), spatial=T),
                    xyFromCell(r, cell=c(S_nodes_cells_r), spatial=T))

r_NS_nodes <- rasterize(x = NS_nodes_r, y = r)

plot(r_NS_nodes)

Except, now after only changing the shape and resolution, the rasterize function returns only the top row! I'm puzzled...

I assume there must be an issue with rasterize, because the r_NS_nodes contains coordinates for both the top and bottom row:

plot(r); points(NS_nodes_r)

Please, if anyone could explain this to me, or provide a solution, I would be greatly appreciative.

-Alex.

Solution

  • The values are there, but not shown on your plot (there are not enough pixels on your screen):

    r_NS_nodes[nrow(r_NS_nodes), 1:3]
#[1] 786 787 788

    Also see

    a <- aggregate(r_NS_nodes, 100, na.rm=T)
plot(a)

    Here is another way to get the cell numbers

    library(raster)
r <- raster(xmn=-150, xmx=-120, ymx=60, ymn=30, ncol=36, nrow=18)
values(r) <- 1:ncell(r)

cells <- cellFromRowColCombine(r, c(1, nrow(r)), 1:ncol(r))

    Now you can do

    xy <- xyFromCell(r, cells)
rr <- rasterize(xy, r)

    or

    rout <- raster(r)
rout[cells] <- 1:length(cells)

    Or you can forego computing cell numbers and do

    rout <- raster(r)
cols <- 1:ncol(r)
rout[c(1, nrow(r)), ] <- c(cols, cols+ncol(r))