0

I have seen a couple of similar questions but they all seem related to situations where numeric values are actually factors. This is not the case here.

I have 2 rasters say A and B. B is cropped from A using raster::rasterize so they have the same resolution, same CRS, etc. I want to find the grid indices where raster B's cells are from raster A. So I use raster::coordinates to pull the coordinate values of both rasters and then try to search for the coordinates of a cell in B within A. The x value seems to work but the y does not yet the values when printed are identical. It only works when I convert them to character.

The code:

p.s. find the rasters here in case you want to test

a.tif = https://github.com/chrisvwn/stackoverflow_coord_match_rasters/blob/master/a.tif

b.tif = https://github.com/chrisvwn/stackoverflow_coord_match_rasters/blob/master/b.tif

> a <- raster::raster("a.tif")
> b <- raster::raster("b.tif")

> a
class      : RasterLayer 
dimensions : 1155, 1441, 1664355  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 2.670833, 14.67917, 4.270834, 13.89583  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 +no_defs 
source     : a.tif 
names      : NL_NGA_OLS.Y_1992_STABLE_LIGHTS.MTSALL.MEAN.RGFF_GADM.3.6.SHPZIP 

> b
class      : RasterLayer 
dimensions : 9, 9, 81  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 7.312499, 7.387499, 5.054167, 5.129167  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 +no_defs 
source     : memory
names      : NL_NGA_OLS.Y_1992_STABLE_LIGHTS.MTSALL.MEAN.RGFF_GADM.3.6.SHPZIP 
values     : 8, 58  (min, max)

> coord_a <- raster::coordinates(a)
> coord_b <- raster::coordinates(b)

> class(coord_a)
"matrix"

> class(coord_b)
"matrix"

> str(coord_a)
 num [1:1664355, 1:2] 2.67 2.68 2.69 2.7 2.71 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:2] "x" "y"

> str(coord_b)
 num [1:81, 1:2] 7.32 7.32 7.33 7.34 7.35 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:2] "x" "y"

> which(coord_a[,1] == coord_b[1,1])
[1]     558    1999    3440    4881    6322    7763    9204   10645   12086   13527   14968   16409   17850
[14]   19291   20732   22173   23614   25055   26496   27937   29378   30819   32260   33701   35142   36583
...
[ reached getOption("max.print") -- omitted 155 entries ]

> which(coord_a[,2] == coord_b[1,2])
integer(0)

> which(as.character(a[,2]) == as.character(b[1,2]))
[1] 1515933 1515934 1515935 1515936 1515937 1515938 1515939 1515940 1515941 1515942 1515943 1515944 1515945
...
[ reached getOption("max.print") -- omitted 441 entries ]

> which(coord_a[,1] == coord_b[1,1] & coord_a[,2] == coord_b[1,2])
integer(0)

> which(coord_a[,1] == coord_b[1,1] & as.character(coord_a[,2]) == as.character(coord_b[1,2]))
[1] 1516490

> coord_a[1516490,]
       x        y 
7.316666 5.125000

> coord_b[1,]
       x        y 
7.316666 5.125000

> print(coord_a[1516490,], digits=10)
          x           y 
7.316665917 5.125000279 

> print(coord_b[1,], digits=10)
          x           y 
7.316665917 5.125000279 

> coord_a[1516490,] == b[1,]
    x     y 
 TRUE FALSE

> as.character(coord_a[1516490,]) == as.character(coord_b[1,])
[1] TRUE TRUE

#for good measure test individually
> coord_a[1516490,1] == coord_b[1,1]
   x 
TRUE 

> coord_a[1516490,2] == coord_b[1,2]
    y 
FALSE 

> as.character(coord_a[1516490,2]) == as.character(coord_b[1,2])
[1] TRUE

> dput(coord_a)
structure(c( ... [very many numeric values]
4.27500028264978, 4.27500028264978, 4.27500028264978, 4.27500028264978, 
4.27500028264978, 4.27500028264978, 4.27500028264978), .Dim = c(1664355L, 
2L), .Dimnames = list(NULL, c("x", "y")))

> dput(coord_b)
structure(c(7.31666591663393, 7.32499924993373, 7.33333258323352, 
7.34166591653331, 7.3499992498331, 7.35833258313289, 7.36666591643269, 
7.37499924973248, 7.38333258303227, 7.31666591663393, 7.32499924993373, 
7.33333258323352, 7.34166591653331, 7.3499992498331, 7.35833258313289, 
7.36666591643269, 7.37499924973248, 7.38333258303227, 7.31666591663393, 
7.32499924993373, 7.33333258323352, 7.34166591653331, 7.3499992498331, 
7.35833258313289, 7.36666591643269, 7.37499924973248, 7.38333258303227, 
7.31666591663393, 7.32499924993373, 7.33333258323352, 7.34166591653331, 
7.3499992498331, 7.35833258313289, 7.36666591643269, 7.37499924973248, 
7.38333258303227, 7.31666591663393, 7.32499924993373, 7.33333258323352, 
7.34166591653331, 7.3499992498331, 7.35833258313289, 7.36666591643269, 
7.37499924973248, 7.38333258303227, 7.31666591663393, 7.32499924993373, 
7.33333258323352, 7.34166591653331, 7.3499992498331, 7.35833258313289, 
7.36666591643269, 7.37499924973248, 7.38333258303227, 7.31666591663393, 
7.32499924993373, 7.33333258323352, 7.34166591653331, 7.3499992498331, 
7.35833258313289, 7.36666591643269, 7.37499924973248, 7.38333258303227, 
7.31666591663393, 7.32499924993373, 7.33333258323352, 7.34166591653331, 
7.3499992498331, 7.35833258313289, 7.36666591643269, 7.37499924973248, 
7.38333258303227, 7.31666591663393, 7.32499924993373, 7.33333258323352, 
7.34166591653331, 7.3499992498331, 7.35833258313289, 7.36666591643269, 
7.37499924973248, 7.38333258303227, 5.12500027920563, 5.12500027920563, 
5.12500027920563, 5.12500027920563, 5.12500027920563, 5.12500027920563, 
5.12500027920563, 5.12500027920563, 5.12500027920563, 5.11666694590606, 
5.11666694590606, 5.11666694590606, 5.11666694590606, 5.11666694590606, 
5.11666694590606, 5.11666694590606, 5.11666694590606, 5.11666694590606, 
5.10833361260649, 5.10833361260649, 5.10833361260649, 5.10833361260649, 
5.10833361260649, 5.10833361260649, 5.10833361260649, 5.10833361260649, 
5.10833361260649, 5.10000027930693, 5.10000027930693, 5.10000027930693, 
5.10000027930693, 5.10000027930693, 5.10000027930693, 5.10000027930693, 
5.10000027930693, 5.10000027930693, 5.09166694600736, 5.09166694600736, 
5.09166694600736, 5.09166694600736, 5.09166694600736, 5.09166694600736, 
5.09166694600736, 5.09166694600736, 5.09166694600736, 5.08333361270779, 
5.08333361270779, 5.08333361270779, 5.08333361270779, 5.08333361270779, 
5.08333361270779, 5.08333361270779, 5.08333361270779, 5.08333361270779, 
5.07500027940822, 5.07500027940822, 5.07500027940822, 5.07500027940822, 
5.07500027940822, 5.07500027940822, 5.07500027940822, 5.07500027940822, 
5.07500027940822, 5.06666694610866, 5.06666694610866, 5.06666694610866, 
5.06666694610866, 5.06666694610866, 5.06666694610866, 5.06666694610866, 
5.06666694610866, 5.06666694610866, 5.05833361280909, 5.05833361280909, 
5.05833361280909, 5.05833361280909, 5.05833361280909, 5.05833361280909, 
5.05833361280909, 5.05833361280909, 5.05833361280909), .Dim = c(81L, 
2L), .Dimnames = list(NULL, c("x", "y")))

I have tried to convert to data.frame with the same outcome

Chris Njuguna
  • 335
  • 1
  • 3
  • 12
  • Try `all.equal`. Since `all.equal` is not vectorized, you can write a helper function. `is.equal <- function(x, y, eps = .Machine$double.eps^0.5) abs(x - y) < eps` – Rui Barradas Aug 22 '19 at 15:22
  • Indeed this is the correct answer. So I guess the question I should be asking is why x works and y doesn't. Meaning maybe there was some very slight recalculation on y when the crop happened. – Chris Njuguna Aug 22 '19 at 16:40
  • Given that the extent of A is a decimal number, you never get exact alignment. Because this question was closed; I cannot give you a full answer, but what you are after might be easier achieved like this: `library(raster); a <- raster(ext=extent(2.670833, 14.67917, 4.270834, 13.89583), nrow=1155, ncol=1441); b <- raster(ext=extent(7.312499, 7.387499, 5.054167, 5.129167), nrow=9, ncol=9); x <- cellFromXY(a, xyFromCell(b, 1:ncell(b)))` – Robert Hijmans Aug 22 '19 at 17:12
  • @RobertHijmans This is brilliant! And much shorter than what I was going to try with the searches! I haven't seen cellFromXY and xyFromCell before so I have to research these. How would I pull out the coordinates as well just so I have exact values that I can use further on? Not sure if this is common knowledge - I am thinking of creating a new question just so I can ask you to put a full answer. – Chris Njuguna Aug 22 '19 at 17:54
  • `xyFromCell` gives you the coordinates. `xyFromCell(b, 1:ncell(b))` is the same as `coordinates(b)`, but here I thought that `xyFromCell` was more instructive. – Robert Hijmans Aug 22 '19 at 18:28
  • Thanks @RobertHijmans makes sense and better than coordinates because the 'xy' functions allow for inverse operations. – Chris Njuguna Aug 22 '19 at 18:41

0 Answers0