Tuesday, 12 January 2016

R Raster: extract weighted mean within circle of specific radius


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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...