2

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?

0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.