2

I am currently using a raster that is a land cover class raster (with 50 different land cover classes at a resolution of 30m). I would like to obtain 50 different rasters (essentially predictors for a future regression analysis) that are data of a continuous nature.

Previously, I used the focal statistics tool in ArcMap but I ended up obtaining either a 0 or a 1 based on areas where a particular land cover type was present.

See image below (All areas of White are 1 and Areas in black are 0):

enter image description here

Is there a way to obtain multiple rasters from a single land cover class raster, with the proportion of each land cover class calculated around a window of 5*5 from a particular point (instead of a 0 or a 1) in R?

I noticed that there is a focal function in the raster package.

  • By nature, the data is nominal. What would a continuous range [0-1] represent? Are you want to calculate the proportion of each landcover class within a specified window? If this is the case then yes, the focal function is what you are after. Please search the site as, I have answered variations of this question at least twice. – Jeffrey Evans Apr 10 '17 at 15:20
  • unless are you intending to change the resolution (coarser), how are you going to change the categorical land cover class to continuous 0-1 values? – Sam Apr 10 '17 at 15:31
  • @JeffreyEvans Thank you for the correction. Yes, I am trying to calculate the proportion of each landcover class within a particular window. – Vijay Ramesh Apr 10 '17 at 15:36
  • @Sam,I am trying to resample it (coarser) to a 1 km to match the extent and resolution of the rest of my predictors. – Vijay Ramesh Apr 10 '17 at 15:36
  • ok thanks, is this 1km res on the same CRS and spatial extent? – Sam Apr 10 '17 at 15:38
4

For calculating proportions around points you can define a function in the raster::extract function.

First, create an example raster with values of [1,2] and generate a random point sample.

library(raster)
library(sp)
 r <- raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, 
             ymx=4453690, resolution=270, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
             +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
 r[] <- rpois(ncell(r), lambda=1)
 r <- calc(r, fun=function(x) { x[x >= 1] <- 1; return(x+1) } )
 x <- sampleRandom(r, 10, na.rm = TRUE, sp = TRUE)
plot(r)
plot(x, add=TRUE, pch=20)

Now we can pull the proportions of each raster class, using a 500m buffer around each point, and store them in a list object. The iterator "i" is defining the raster class value in the loop and the function being passed to raster::extract.

 landcover.prop <- list()
   for( i in 1:2) {
     landcover.prop[[i]] <- extract(r, x, buffer=500, small=TRUE, fun=function(x, p = i) 
                                   { prop.table( ifelse(x == p, 1, 0) ) } )
    }

The elements in the list object are vectors that are ordered the same as your points so, you can just add them to the SpatialPointsDataFrame object.

  • Thank you @Jeffrey Evans. This is precisely what I was looking for. To be clear, the iterator i is moving through the different land cover classes as you mentioned right? Currently running your example code. – Vijay Ramesh Apr 11 '17 at 12:37
1

based on the info you have given, i hope this helps (R):

# your 30m raster with various land cover classes (here dummy data with 5)
r <- raster(ncol=10,nrow=10)
r[] <- sample(seq(from = 1, to = 5, by = 1), size = 100, replace = TRUE)

# create an empty raster stack ready to receive rasters
st <- stack() 

# Then, loop through the land cover classes, convert to binary (at original resolution) 
# and aggregate up to new desired resolution.

for(i in 1:maxValue(r)){
  d <- r
  d[d!=i] <- NA
  d[!is.na(d)] <- 1
  dagg <- aggregate(d,fact=c(5,5),fun=sum)
  print(cellStats(dagg,sum))
}

The result of each cellStats(sum) should be the same as freq(r).

However, if you are indeed doing 1km aggregation or 'window' over 30m data, which obviously isn't a factor, then you will have to consider converting the 1km summary raster into polygons and doing some summary stats instead (more involved, a separate answer i think).

  • Thanks, this is definitely an interesting solution. In essence, this would give me the total proportion of land-cover at every 5 *5 window right? I think I might post a separate question again, as I have begun to think about the same, but would love to obtain the proportion of each land cover class around a given point, in a 5 * 5 km window. – Vijay Ramesh Apr 10 '17 at 16:02
  • Hi @Sam, I have modified the question to essentially get at the proportion of each land cover class around a given point. – Vijay Ramesh Apr 10 '17 at 20:24
  • @Vijay, i can't answer now, but essentially what my very simple example did was subset the raster for 1 land class, turn it binary (1 = land class exists, NA it doesn't), then aggregated to a desired new resolution (my example was 25 times coarser). Then you can divide the new raster by the aggregation factor to get a 0 to 1 value (i didnt do this above). – Sam Apr 10 '17 at 20:33
  • @Vijay, you have a slightly different qu now though, a moving focal window? is this what you want? to me that sounds very odd but if it's what you want.... – Sam Apr 10 '17 at 20:34
  • 3
    Here is where you really need to be much clearer. In R the method for extracting data around a point is much different than creating new rasters representing proportion(s). Please look at ?extract. You can pass a custom function and specify a window around each point. This is computationally, much more efficient than deriving new rasters using focal. – Jeffrey Evans Apr 10 '17 at 20:37

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.