Collectives™ on Stack Overflow

Find centralized, trusted content and collaborate around the technologies you use most.

Teams

Q&A for work

Connect and share knowledge within a single location that is structured and easy to search.

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.
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) 
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))
        
Thanks for contributing an answer to Stack Overflow!
Please be sure to answer the question. Provide details and share your research!
But avoid …
Asking for help, clarification, or responding to other answers.
Making statements based on opinion; back them up with references or personal experience.
To learn more, see our tips on writing great answers.