0

I have a raster of 1 km wind direction data that I need to run a focal buffer on. Instead of using the default buffering mean function, I need to use vector averaging since my wind direction is in degrees (e.g. (360+0)/2 = 180/South, instead of 360/North!).

I tried using the circular mean function from the circular package, but this doesn't work with rasters.

How do I perform vector average focal buffers on rasters?

I wrapped avgWind function around circular mean but still get the same error.

library(raster)
library(circular)

# creating a sample raster for stackoverflow
xy <- matrix(sample(0:359), 10, 10)
r <- raster(xy)

avgWind <- function(r, units) {
  r_avg <- mean(circular(r, units = "degrees"))
  return(r_avg)
}

r_buf <- focal(r, focalWeight(r, d=20, type='circle')
               ,fun = avgWind(r)
               ,na.rm = T)


Error in x/180 : non-numeric argument to binary operator
In addition: Warning message:
In class(x) <- c("circular", cl) :
  Setting class(x) to multiple strings ("circular", "RasterLayer", ...); result will no longer be an S4 object

EDIT

After applying obrl_soil's example to my projected data the buffered data appears to be inverted.

10k buffered wind direction raster

original wind direction raster

  • 1
    Your fun argument needs to be an actual function. – Jeffrey Evans Nov 28 '18 at 4:18
  • I still get the same error after wrapping circular into its own function. – philiporlando Nov 28 '18 at 4:25
  • 1
    I'd usually decompose to horizontal and vertical components of the wind and average those, then calculate the direction again. Do you have magnitude as well as direction? – mdsumner Nov 28 '18 at 5:17
  • 1
    The syntax problem needed to fix here is delete "(r)" from "avgWind(r)" - it wants the function, not a return value from it. (But running this code crashes R for me so I'd check what avgWind is returning) – mdsumner Nov 28 '18 at 5:22
  • My 32GB workstation could handle it, but it's crashing on my 4GB laptop too. I'll look into calculating horizontal and vertical components of wind direction/magnitude and see what I can do. – philiporlando Nov 28 '18 at 5:50
2

focalWeight appears designed to work on a projected raster, so I think half the problems you're having actually come from the reprex you've built. Terrible irony! To build a circular filter for unprojected data, I like the method in this blog post:

make_circ_filter <- function(radius, res) {
  circ_filter<-matrix(0, nrow=1+(2*radius/res), 
                         ncol=1+(2*radius/res))
  dimnames(circ_filter)[[1]]<-seq(-radius, radius, by=res)
  dimnames(circ_filter)[[2]]<-seq(-radius, radius, by=res)
  sweeper<-function(mat){
    for(row in 1:nrow(mat)){
      for(col in 1:ncol(mat)){
        dist<-sqrt((as.numeric(dimnames(mat)[[1]])[row])^2 +
                     (as.numeric(dimnames(mat)[[1]])[col])^2)
        if(dist<=radius) {mat[row, col]<-1}
      }
    }
    return(mat)
  }
  out<-sweeper(circ_filter)
  return(out)
}


cf <- make_circ_filter(5, 1)
# for an efficient mean filter, adjust the non-NA values 
# in cf to 1/n where n is the number of non-NA cells, and
# then set NA to 0
cf[which(cf == 1)] <- 1/length(cf[!is.na(cf)])
cf[which(is.na(cf))] <- 0

EDIT

Ok the technical note you've supplied has a much better solution. So here's how I would spatialise it:

library(raster)
library(slga)

# handy data
aoi <-  c(148.3, -21.2, 148.7, -20.8)
aspect <- slga::get_lscape_data('ASPCT', aoi )
fake_speed <- slga::get_lscape_data('SLPPC', aoi) # its really slope >.>

# then follow tech note method to average wind
rad2deg <- function(rad) {(rad * 180) / (pi)}
deg2rad <- function(deg) {(deg * pi) / (180)}

# calculate components u and v
dir_comp_u <- fake_speed * sin(deg2rad(aspect))
dir_comp_v <- fake_speed * cos(deg2rad(aspect))

# run focal on each component to get a spatial mean
avg_comp_u <- focal(dir_comp_u, cf)
avg_comp_v <- focal(dir_comp_v, cf)

# back-convert
wd_avg_rad <- calc(stack(avg_comp_u, avg_comp_v), 
                   function(cell) {
                     if(any(is.na(cell))) { 
                       NA_real_
                     } else {
                         atan2(cell[1], cell[2])
                     }})

# and then to degrees
wd_avg_deg <- abs(rad2deg(wd_avg_rad) - 180)

I don't really know how to interpret the output as I'm using fake data, but it seems reasonable...?

  • Wow this is awesome! Thank you! Would you recommend adding in the wind speed component and modifying your code to follow this example? researchgate.net/profile/Stuart_Grange2/publication/… – philiporlando Nov 28 '18 at 8:08
  • 1
    And I'd use this method to spatially smooth aspect data as well - just leave the speed component out. – obrl_soil Nov 28 '18 at 9:53
  • The spatial buffers are working, but the final output appears to be inverted. I've added comparison figures of my projected data, swapping cf for a focalWeight parameter. I've gone through the math several times and haven't figured out where the problem is yet. – philiporlando Nov 29 '18 at 1:20
  • 1
    Ditch the 180 in the last step, I guess? Exactly what to do depends on whether your origin point is at 12 o'clock or elsewhere. Different GIS apps set different origins. – obrl_soil Nov 29 '18 at 4:19
  • 1
    I was missing a negative sign when calculating dir_comp_u and dir_comp_v. It only took me a few weeks to figure it out ha! Thanks again for all of your help on this! – philiporlando Dec 24 '18 at 6:10

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.