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

为什么光栅提取会产生暗淡的误差?

  •  0
  • Lukas  · 技术社区  · 6 年前

    我试图从光栅砖中提取值并得到平均值,但得到的错误我认为与光栅砖的尺寸有关。

    NOAA

    library(raster)
    
    ERSST <- rotate(brick('sst.mnmean.nc'))
    ERSST
    class       : RasterBrick 
    dimensions  : 89, 180, 16020, 1976  (nrow, ncol, ncell, nlayers)
    resolution  : 2, 2  (x, y)
    extent      : -179, 181, -89, 89  (xmin, xmax, ymin, ymax) # ignore extent, needs fixing but not relevant for the question
    coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
    data source : in memory
    names       : X1854.01.01, X1854.02.01, X1854.03.01, X1854.04.01, X1854.05.01, X1854.06.01, X1854.07.01, X1854.08.01, X1854.09.01, X1854.10.01, X1854.11.01, X1854.12.01, X1855.01.01, X1855.02.01, X1855.03.01, ... 
    min values  :        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8,        -1.8, ... 
    max values  :    32.09937,    31.44189,    31.72137,    31.47466,    33.23633,    31.90788,    35.22922,    33.66898,    32.26702,    32.15502,    31.68270,    31.74512,    31.32458,    31.23049,    29.88974, ... 
    Date        : 1854-01-01, 2018-08-01 (min, max)
    

    任意点

    xy <- data.frame(x = -49, y = 45)
    

    extract(ERSST, xy, buffer = 1e+05, small = TRUE, fun = mean)
    Error in apply(x, 2, fun2) : dim(X) must have a positive length
    

    我之所以认为这是一个维度问题,是因为我在尝试指定要使用的层时遇到了错误

    extract(ERSST, xy, buffer = 1e+05, small = TRUE, layer = 10, nl = 10)
    Error in x[, lyrs] : incorrect number of dimensions
    

    mERSST <- mean(ERSST)
    extract(mERSST, xy, buffer = 1e+05, small = TRUE, fun = mean)
    [1] 5.649212
    

    @RobertHijmans的答案让我意识到,我总是从extract中只得到一个值,即使点位于几个网格单元的连接处,如上面的示例所示。

    plot(MSST, xlim = c(-60, -40), ylim = c(40, 50))
    points(xy)
    

    enter image description here

    extract(mERSST, xy, buffer = 1e+05, small = TRUE, cellnumbers = TRUE)
    [[1]]
           cell       value 
    3845.000000    5.649212 
    

    我只得到一个值,而我希望有4个,不管缓冲区有多小。我是不是遗漏了什么?所以我试着把我的点转换成一个圆,然后用它来提取数据

    coordinates(xy) <- ~ x + y
    proj4string(xy) <- '+init=epsg:4326'
    
    xy_utm <- spTransform(xy, CRS('+init=epsg:32621'))
    gbf_utm <- rgeos::gBuffer(xy_utm, width = 1e5, quadsegs = 250L)
    gbf <- spTransform(gbf_utm, CRS(proj4string(xy)))
    
    
    plot(ERSST[[1]], xlim = c(-60, -40), ylim = c(40, 50))
    points(xy, pch = 19)
    plot(gbf, add = TRUE)
    

    enter image description here

    extract(ERSST[[1]], gbf, small = TRUE, weights = TRUE)
    

    [[1]]
            value weight
    [1,] 1.722664   0.25
    [2,] 3.683457   0.25
    [3,] 5.985203   0.25
    [4,] 8.442450   0.25
    

    在版本2.6.7中(这似乎是有道理的)。

    但是

    [[1]]
              value      weight
      [1,] 1.722664 0.001236928
      [2,] 1.722664 0.003935680
      [3,] 1.722664 0.005285056
      [4,] 3.683457 0.005285056
      [5,] 3.683457 0.003935680
      [6,] 3.683457 0.001236928
      [7,] 1.722664 0.002136512
      [8,] 1.722664 0.008321151
      [9,] 1.722664 0.011244799
     [10,] 1.722664 0.011244799
     [11,] 1.722664 0.011244799
     [12,] 3.683457 0.011244799
     [13,] 3.683457 0.011244799
     [14,] 3.683457 0.011244799
     [15,] 3.683457 0.008208703
     [16,] 3.683457 0.001911616
     [17,] 1.722664 0.003036096
     [18,] 1.722664 0.010907455
     [19,] 1.722664 0.011244799
     [20,] 1.722664 0.011244799
     [21,] 1.722664 0.011244799
     [22,] 1.722664 0.011244799
     [23,] 3.683457 0.011244799
     [24,] 3.683457 0.011244799
     [25,] 3.683457 0.011244799
     [26,] 3.683457 0.011244799
     [27,] 3.683457 0.010907455
     [28,] 3.683457 0.003036096
     [29,] 1.722664 0.000449792
     [30,] 1.722664 0.010232767
     [31,] 1.722664 0.011244799
     [32,] 1.722664 0.011244799
     [33,] 1.722664 0.011244799
     [34,] 1.722664 0.011244799
     [35,] 1.722664 0.011244799
     [36,] 3.683457 0.011244799
     [37,] 3.683457 0.011244799
     [38,] 3.683457 0.011244799
     [39,] 3.683457 0.011244799
     [40,] 3.683457 0.011244799
     [41,] 3.683457 0.010232767
     [42,] 3.683457 0.000337344
     [43,] 1.722664 0.003036096
     [44,] 1.722664 0.011244799
     [45,] 1.722664 0.011244799
     [46,] 1.722664 0.011244799
     [47,] 1.722664 0.011244799
     [48,] 1.722664 0.011244799
     [49,] 1.722664 0.011244799
     [50,] 3.683457 0.011244799
     [51,] 3.683457 0.011244799
     [52,] 3.683457 0.011244799
     [53,] 3.683457 0.011244799
     [54,] 3.683457 0.011244799
     [55,] 3.683457 0.011244799
     [56,] 3.683457 0.002923648
     [57,] 5.985203 0.002923648
     [58,] 5.985203 0.011244799
     [59,] 5.985203 0.011244799
     [60,] 5.985203 0.011244799
     [61,] 5.985203 0.011244799
     [62,] 5.985203 0.011244799
     [63,] 5.985203 0.011244799
     [64,] 8.442450 0.011244799
     [65,] 8.442450 0.011244799
     [66,] 8.442450 0.011244799
     [67,] 8.442450 0.011244799
     [68,] 8.442450 0.011244799
     [69,] 8.442450 0.011244799
     [70,] 8.442450 0.002923648
     [71,] 5.985203 0.000337344
     [72,] 5.985203 0.010120319
     [73,] 5.985203 0.011244799
     [74,] 5.985203 0.011244799
     [75,] 5.985203 0.011244799
     [76,] 5.985203 0.011244799
     [77,] 5.985203 0.011244799
     [78,] 8.442450 0.011244799
     [79,] 8.442450 0.011244799
     [80,] 8.442450 0.011244799
     [81,] 8.442450 0.011244799
     [82,] 8.442450 0.011244799
     [83,] 8.442450 0.010007871
     [84,] 8.442450 0.000224896
     [85,] 5.985203 0.002811200
     [86,] 5.985203 0.010795007
     [87,] 5.985203 0.011244799
     [88,] 5.985203 0.011244799
     [89,] 5.985203 0.011244799
     [90,] 5.985203 0.011244799
     [91,] 8.442450 0.011244799
     [92,] 8.442450 0.011244799
     [93,] 8.442450 0.011244799
     [94,] 8.442450 0.011244799
     [95,] 8.442450 0.010682559
     [96,] 8.442450 0.002698752
     [97,] 5.985203 0.001799168
     [98,] 5.985203 0.007871359
     [99,] 5.985203 0.011244799
    [100,] 5.985203 0.011244799
    [101,] 5.985203 0.011244799
    [102,] 8.442450 0.011244799
    [103,] 8.442450 0.011244799
    [104,] 8.442450 0.011244799
    [105,] 8.442450 0.007871359
    [106,] 8.442450 0.001799168
    [107,] 5.985203 0.001236928
    [108,] 5.985203 0.003935680
    [109,] 5.985203 0.005285056
    [110,] 8.442450 0.005285056
    [111,] 8.442450 0.003935680
    [112,] 8.442450 0.001236928
    

    在版本2.7-13中,这不可能是正确的。

    2 回复  |  直到 6 年前
        1
  •  1
  •   Robert Hijmans    6 年前

    为了用一种简化的方式来说明上面的讨论(在解决了多边形提取问题之后),我得到

    library(raster)
    
    r <- raster(xmn=-59, xmx=-39, ymn=41, ymx=49, res=2, vals=1:40)
    
    xy <- SpatialPoints(data.frame(x = -49, y = 45), proj4string = CRS('+init=epsg:4326'))
    p <- buffer(xy, width = 1e5, quadsegs = 250L)
    plot(r)
    plot(p, add=T)
    
    extract(r, xy)       
    #26 
    extract(r, p)
    #[[1]]
    #[1] 16 26 25 15
    
    extract(r, p, weights=T)
    #[[1]]
    #     value weight
    #[1,]    15   0.25
    #[2,]    16   0.25
    #[3,]    25   0.25
    #[4,]    26   0.25
    
    extract(r, xy, buffer=100000)
    #[[1]]
    #value 
    #   15 
    
    extract(r, xy, buffer=1000000)
    #[[1]]
    # [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
    #[37] 37 38 39 40
    
        2
  •  1
  •   www    6 年前

    layer nl

    一种解决方法是先提取值,然后对值进行子集。

    library(raster)
    
    value <- extract(ERSST, point, buffer = 1e+05, small = TRUE)
    value[[1]][10:19]
    # X1854.10.01 X1854.11.01 X1854.12.01 X1855.01.01 X1855.02.01 X1855.03.01 X1855.04.01 X1855.05.01 
    # 7.932416    5.712043    4.428292    3.010927    2.289096    2.385752    2.528488    3.261783 
    # X1855.06.01 X1855.07.01 
    # 5.762860    7.617740