I'm trying to find the index of a specific location in a spatial raster. But depending on how I calculate the index I'm getting different answers, and I can't figure out why they're different or which one is correct.
Please find below, a minimal working example:
> library(raster)
> library(tidyverse)
>
> ### define raster (grid) & point of interest (poi)
> lon <- seq(-124.6875, -114.0625, by = 0.625)
> lat <- seq(32.5, 49.0, by = 0.5)
> grid <- expand.grid(lon,lat) %>%
+ mutate(values = rnorm(nrow(.))) %>%
+ rasterFromXYZ(crs = "+proj=longlat +datum=WGS84 +no_defs")
> poi <- c(-122.8125, 38.5)
>
> ### option 1: find row & column indices manually ###
>
> # row-col index:
> c <- which(lon==poi[1]) #longitude = x = matrix columns
> r <- which(lat==poi[2]) #latitude = y = matrix rows
> c(row=r, col=c)
row col
13 4
>
> # cell index:
> cellFromRowCol(grid, row=r, col=c)
[1] 220
>
> # cell value:
> grid[r,c]
[1] -0.6369638
>
> ### option 2: use raster functions ###
>
> # cell index:
> cell <- cellFromXY(grid, xy = poi)
> cell
[1] 382
>
> # row-col index:
> rowColFromCell(grid, cell)
row col
[1,] 22 4
>
> # cell value:
> grid[cell]
[1] -0.8445352
>
> ## none of these match the results from the manual method!
As it stands, Option 1 is incorrect and Option 2 is correct.
The error in Option 1 comes from the fact that, in a raster, cells are numbered from the upper left cell to the upper right cell and then continuing on the left side of the next row, and so on until the last cell at the lower-right side of the raster. So, if you want to use Option 1 and get the right result (i.e. the same as Option 2) you need to reverse the values of your lat vector (i.e. set the values in decreasing order).
Please find below a reprex.
Reprex
library(raster)
library(tidyverse)
### define raster (grid) & point of interest (poi)
lon <- seq(-124.6875, -114.0625, by = 0.625)
lat <- rev(seq(32.5, 49.0, by = 0.5)) # !!!! Reverse the order of the 'lat' vector
set.seed(1234) # !!!!Set the seed for the reprex to be reproducible
grid <- expand.grid(lon,lat) %>%
mutate(values = rnorm(nrow(.))) %>%
rasterFromXYZ(crs = "+proj=longlat +datum=WGS84 +no_defs")
poi <- c(-122.8125, 38.5)
### option 1: find row & column indices manually ###
# row-col index:
c <- which(lon==poi[1]) #longitude = x = matrix columns
r <- which(lat==poi[2]) #latitude = y = matrix rows
c(row=r, col=c)
#> row col
#> 22 4
# cell index:
cellFromRowCol(grid, row=r, col=c)
#> [1] 382
# cell value:
grid[r,c]
#> [1] -2.651741
### option 2: use raster functions ###
# cell index:
cell <- cellFromXY(grid, xy = poi)
cell
#> [1] 382
# cell value
grid[cell]
#> [1] -2.651741
Created on 2022-03-24 by the reprex package (v2.0.1)