I have a raster object from which I want to extract area weighted mean values from a circle with a specific radius in m. The raster has latitude and longitude wrt to the WGS84 datum.
I can extract values within a circle using extract
and specifying a buffer
, but how can I get a weighted mean?
The raster layer can be downloaded here.
WOA <- readRDS('WOA.RDS')
crs(WOA)
CRS arguments:
+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
xy <- data.frame(x = -40, y = 60)
# extract values from cells from circle with 100 km radius
extract(WOA, xy, buffer = 1e+5)
but to make use of the weights
in extract
I think need to convert xy to a SpatialPoint
xy_sp <- SpatialPoints(xy, proj4string = CRS('+proj=longlat +datum=WGS84'))
and use gBuffer
from the rgeos
package to expand the point to a circle with the desired radius (i.e. a SpatialPolygon), but how do I take care of the units since xt_sp
is in degrees? Can I use spTransform
(rgdal) and what (metric) coordinate system would be appropriate for my global dataset?
Answer
I usually spTransform()
the point data, for which I need a buffer in meters, to UTM which is based on +units=m
. Yours falls into UTM zone 21N (ie. EPSG:32621); if you're uncertain about the right UTM zone, feel free to use online tools like this that easily provide you with the required information.
Then, create a 100-km buffer from rgeos::gBuffer()
and spTransform()
the resulting 'SpatialPolygons' object back to Lat/Long (ie. EPSG:4326) in order to use it as y-input for extract()
along with your 'Raster*' object.
## transform point data to 'SpatialPoints'
coordinates(xy) <- ~ x + y
proj4string(xy) <- "+init=epsg:4326"
## create 100-km buffer
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)))
## extract values and calculate weighted mean
vls = extract(WOA, gbf, weights = TRUE)[[1]]
sum(vls[, 1] * vls[, 2])
# [1] 4.18117
No comments:
Post a Comment