Reputation: 727
I'm trying to extract values and get the mean from a raster brick, but get an error I think is related to the dimensions of the raster brick.
the data have been downloaded from NOAA
what I've done is the following:
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)
Arbitrary point
xy <- data.frame(x = -49, y = 45)
and when I extract I get:
extract(ERSST, xy, buffer = 1e+05, small = TRUE, fun = mean)
Error in apply(x, 2, fun2) : dim(X) must have a positive length
what makes me think it's a dimension problem is the error I get when I try to specify the layers to use
extract(ERSST, xy, buffer = 1e+05, small = TRUE, layer = 10, nl = 10)
Error in x[, lyrs] : incorrect number of dimensions
and it seems to work fine if I average first (but this is not what I want, I need the time series at the point)
mERSST <- mean(ERSST)
extract(mERSST, xy, buffer = 1e+05, small = TRUE, fun = mean)
[1] 5.649212
Perhaps it's the Date attribute in the raster brick. Any workarounds or solutions to prevent this error?
The answer by @RobertHijmans made me realise that I always get just a single value from extract, even when the point is at the junction of several grid cells, as in the example above.
plot(MSST, xlim = c(-60, -40), ylim = c(40, 50))
points(xy)
Using:
extract(mERSST, xy, buffer = 1e+05, small = TRUE, cellnumbers = TRUE)
[[1]]
cell value
3845.000000 5.649212
I get only a single value, whereas I would expect there to be 4, no matter how small the buffer. Am I missing something in extract? So I tried with converting my point to a circle and use that to extract data
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)
extract(ERSST[[1]], gbf, small = TRUE, weights = TRUE)
this gives me:
[[1]]
value weight
[1,] 1.722664 0.25
[2,] 3.683457 0.25
[3,] 5.985203 0.25
[4,] 8.442450 0.25
in version 2.6.7 (and this seems to make sense).
but
[[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
in version 2.7-13, which can't be right.
Upvotes: 0
Views: 600
Reputation: 47536
To illustrate the discussion above in a simplified manner (and after fixing the extract with polygons issue), I get
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
Upvotes: 1
Reputation: 39174
It looks like the layer
and nl
arguments are not working. I am not sure why.
One workaround is to extract values first and then subset the values.
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
Upvotes: 1