1

I want to find the values of row and column present in the cell. I have written a code to fetch value as given below:

x<-raster()
    values(x)<-1:ncell(x)
class       : RasterLayer 
dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 1, 64800  (min, max)

There is another file as given below which contains lat and long information. So, we have to find these points in the raster object and then store values in the format given at last.

long<-c(106.61291,-81.97224,-84.4277,-97.66631,-72.68604,-73.1247,-119.05662,
   -86.75413,-108.5377,-86.67737,-116.22231,-71.00956,-91.15146,-73.1516)
lat<-c(35.04333,33.37378,33.64073,30.19743,41.93887,41.1642,35.43291,33.56243,
       45.8037,36.12632,43.56582,42.36561,30.53236,44.4707)

xy<-data.frame(long,lat)
coordinates(xy)=c('long','lat')
proj4string(xy)=proj4string(x)
points<-cellFromXY(x,xy)

so , points contains all the cell number where xy is present in raster object x.

As given below:
[1] 19727 20259 20256 21323 17388 17387 19501 20254 15912 19174 16624 17029 21329 16307

I want to find the value of row and column for every cell given in points as I have done below: For example taking cell number 19727 from "points"

rowColFromCell(x,19727)
  row col
  [1,]  55 287

  getValues(x,row=55)[287]
  value:16109

And I want values of 9*9 grid where the given row and column i.e 55,287 is in centre.

And I want to do the same for all the cells given in the points.

Final output should look like this:

 long        lat    cell_1  cell_2  cell_3 cell _4 ......... cell_9
106.61291 35.04333     16109  ...     ....   ....              ....

Thanks in advance!

User0590
  • 103
  • 12
  • I think `getValues(r,row=55)[287]` is 19727. And in the beginning, you created the raster layer as `x`, but at the end, it becomes `r`. – www Sep 10 '17 at 14:56

1 Answers1

3

We can define a function to get the coordinates of the 9*9 grid for each point based on the resolution of the raster layer. After that, we can use the extract function from the raster package to extract the values. In the following example, xy2 is the final output. Column 5 is the value from the center cell. All the other columns show the value from the adjacent 9*9 grid.

library(raster)

# Create the example raster
r <- raster()
values(r) <- 1:ncell(r)

# Create a data frame with long-lat coordinates
long<-c(106.61291,-81.97224,-84.4277,-97.66631,-72.68604,-73.1247,-119.05662,
        -86.75413,-108.5377,-86.67737,-116.22231,-71.00956,-91.15146,-73.1516)
lat<-c(35.04333,33.37378,33.64073,30.19743,41.93887,41.1642,35.43291,33.56243,
       45.8037,36.12632,43.56582,42.36561,30.53236,44.4707)

xy <- data.frame(long, lat)

# Design a function to calculate the coordinates of the surrounding eight cells
coord_nine <- function(long, lat, res){
  cell1 <- c(long - res[1], lat + res[2])
  cell2 <- c(long, lat + res[2])
  cell3 <- c(long + res[1], lat + res[2])
  cell4 <- c(long - res[1], lat)
  cell5 <- c(long, lat)
  cell6 <- c(long + res[1], lat)
  cell7 <- c(long - res[1], lat - res[2])
  cell8 <- c(long, lat - res[2])
  cell9 <- c(long + res[1], lat - res[2])

  coord_list <- list(cell1, cell2, cell3, cell4, cell5, 
                     cell6, cell7, cell8, cell9)

  coord <- do.call(rbind, coord_list)
  return(coord)
}

# Use lapply to loop through all points to get the 9-cell coordinates for each of them
coord_list <- lapply(1:nrow(xy), function(i){
  coord_nine(xy[i, 1], xy[i, 2], res = res(r))
})

# Extract the values based on the coordinates in coord_list
ex_values <- sapply(coord_list, function(coords){
  extract(r, coords)
})

# Combine the results to the original xy data frame
xy2 <- cbind(xy, t(ex_values))

# View the results
xy2
         long      lat     1     2     3     4     5     6     7     8     9
1   106.61291 35.04333 19366 19367 19368 19726 19727 19728 20086 20087 20088
2   -81.97224 33.37378 19898 19899 19900 20258 20259 20260 20618 20619 20620
3   -84.42770 33.64073 19895 19896 19897 20255 20256 20257 20615 20616 20617
4   -97.66631 30.19743 20962 20963 20964 21322 21323 21324 21682 21683 21684
5   -72.68604 41.93887 17027 17028 17029 17387 17388 17389 17747 17748 17749
6   -73.12470 41.16420 17026 17027 17028 17386 17387 17388 17746 17747 17748
7  -119.05662 35.43291 19140 19141 19142 19500 19501 19502 19860 19861 19862
8   -86.75413 33.56243 19893 19894 19895 20253 20254 20255 20613 20614 20615
9  -108.53770 45.80370 15551 15552 15553 15911 15912 15913 16271 16272 16273
10  -86.67737 36.12632 18813 18814 18815 19173 19174 19175 19533 19534 19535
11 -116.22231 43.56582 16263 16264 16265 16623 16624 16625 16983 16984 16985
12  -71.00956 42.36561 16668 16669 16670 17028 17029 17030 17388 17389 17390
13  -91.15146 30.53236 20968 20969 20970 21328 21329 21330 21688 21689 21690
14  -73.15160 44.47070 15946 15947 15948 16306 16307 16308 16666 16667 16668
www
  • 38,575
  • 12
  • 48
  • 84
  • Thanks @ycw, I am getting below error when I am extracting values: `Error in UseMethod("extract_") : no applicable method for 'extract_' applied to an object of class "c('RasterLayer', 'Raster', 'BasicRaster')"`. Can you help me with this? Thanks once again for your help. – User0590 Sep 11 '17 at 01:47
  • @User0590 Did you get the error for my example code or for your actual data? – www Sep 11 '17 at 01:49
  • @User0590 This post could be relevant. https://stackoverflow.com/questions/36616254/rscript-why-is-error-in-usemethodextract-being-indicated-when-attemptin – www Sep 11 '17 at 01:50
  • one more question, in my case `xy` is not a data frame, it's a class of SpatialPoints. So, when I am doing `nrow(xy)`, it's giving `NULL`. Any comments on this? – User0590 Sep 11 '17 at 02:07
  • 1
    @User0590 Based on the post I found, it is likely that you loaded other packages with another `extract` function. You may want to restart your R and make sure only load the `raster` package. Or you may want to add `raster::extract` to make sure you are using the `extract` function from the `raster` package. – www Sep 11 '17 at 02:26
  • @User0590 You may also want to convert your spatial point object to a data frame. – www Sep 11 '17 at 02:28
  • @User0590 You can use `as.data.frame(xy)` to do that. – www Sep 11 '17 at 02:31
  • if I am converting spatial point class to data frame, extract value is Null. – User0590 Sep 11 '17 at 10:42
  • @User0590 Not sure what do you mean. Please make sure the code works for the example I provided. And then figure out what is the difference between my example and your real world data. – www Sep 11 '17 at 10:43
  • Thanks @ycw for your patience.. In my original data, xy looks like this `class : SpatialPoints features : 100 extent : -123.0028, -70.3097, 25.79586, 47.45025 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0`. So, when I am running 'extract(r,coords)', it is giving Null value. I don't know why this is happening. Any comments on this. – User0590 Sep 11 '17 at 13:39
  • And if I am converting xy into data frame, then because of coord ref is not same , extract function is not able to find xy coordinates in raster object r. Because of which it is giving null. – User0590 Sep 11 '17 at 13:42
  • @User0590 If your raster is in a different coordinate system, you should convert your raster's coordinate system to be the same as the `xy`. Or you should convert your `xy` to have the same coordinate system as the raster. – www Sep 11 '17 at 13:51
  • I converted xy to the same coordinate as r coordinates(xy)=c('long','lat') proj4string(xy)=proj4string(r). But then coord_list <- lapply(1:nrow(xy), function(i){ coord_nine(xy[i, 1], xy[i, 2], res = res(r)) }), I am getting error in nrow(xy), xy[i, 1] and xy[i, 2]. – User0590 Sep 11 '17 at 14:21
  • @User0590 I think without access to your actual point data and raster layer, it is difficult to see what happened. I suggest that, if the code works for the example dataset in this post, and if you think my solution is helpful and constructive, accept or upvote my answer. After that, start a new post with reproducible example datasets of your real world data. No need to repeat the all the details in this question here, but describing why the `extract` function from the `raster` is not working. By doing this, you will have a higher chance to get the help you need. – www Sep 11 '17 at 14:43