代码之家  ›  专栏  ›  技术社区  ›  User0590

使用R中的光栅包查找单元格中存在的相邻行和列的值

  •  1
  • User0590  · 技术社区  · 7 年前

    我想找到单元格中存在的行和列的值。我编写了一个代码来获取如下所示的值:

    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)
    

    下面给出的另一个文件包含lat和long信息。因此,我们必须在光栅对象中找到这些点,然后以给定的格式存储值。

    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)
    

    因此,点包含xy存在于光栅对象x中的所有单元号。

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

    例如,从“点”中提取单元格编号19727

    rowColFromCell(x,19727)
      row col
      [1,]  55 287
    
      getValues(x,row=55)[287]
      value:16109
    

    我想要9*9网格的值,其中给定的行和列,即55287在中心。

    我想对点中给出的所有单元格做同样的处理。

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

    提前感谢!

    1 回复  |  直到 7 年前
        1
  •  3
  •   www    7 年前

    我们可以定义一个函数,根据光栅层的分辨率获得每个点的9*9网格坐标。之后,我们可以使用 extract 来自的函数 raster 包以提取值。在以下示例中, xy2 是最终输出。柱 5 是来自中心单元格的值。所有其他列显示来自相邻9*9网格的值。

    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