7

I've been running into all sorts of issues using ArcGIS ZonalStats and thought R could be a great way. Saying that I'm fairly new to R, but got a coding background.

The situation is that I have several rasters and a polygon shapefile with many features of different sizes (though all features are bigger than a raster cell and the polygon features are aligned to the raster). I've figured out how to get the mean value for each polygon feature using the raster library with extract:

#load packages required
require(rgdal)
require(sp)
require(raster)
# ---Set the working directory-------
datdir <- "/test_data/"

#Read in grid of water depth
ras <- raster("test_data/raster/pl_sm_rp1000/w001001.adf")

#read in polygon shape file
proxNA <- shapefile("test_data/proxy/PL_proxy_WD_NA_test") 

#calc mean depth per polygon feature
#unweighted - only assigns grid to district if centroid is in that district
proxNA$RP1000 <- extract(ras, proxNA, fun = mean, na.rm = TRUE, weights = FALSE)

#plot depth values 
spplot(proxNA[,'RP1000'])

The issue I have is that I also need an area based ratio between the area of the polygon and all non NA cells in the same polygon. I know what the cell size of the raster is and I can get the area for each polygon, but the missing link is the count of all non-NA cells in each feature. I managed to get the cell number of all the cells in the polygon proxNA@data$Cnumb1000 <- cellFromPolygon(ras, proxNA)and I'm sure there is a way to get the actual value of the raster cell, which then requires a loop to get the number of all non-NA cells combined with a count, etc. BUT, I'm sure there is a much better and quicker way to do that! If any of you has an idea or can point me in the right direction, I would be very grateful!

  • If I was really done w/debugging the zonalstats approach (which would likely be the ideal way), I would look at numpy before R. That said, is ras holding NA flags using a legit value? Seems like you could either filter for that value or get a count of those values after the fact. – Roland Jan 27 '14 at 17:28
  • @Roland: Thanks! NA is NA in ras and hasn't got a specific value. So what you are saying is that I could filter for NA (or replacement value) and categorise for each polygon to get the count, then subtract from overall cell number. Interesting, but a bit still a bit long-winded. I was hoping for a Count function or somthing along that line. – Hubert Jan 27 '14 at 17:40
  • You could get a seed bump by not using the raster format. The advantage to raster is that it is memory safe. Since you are creating an sp object then coercing to raster you lose the advantage. Keeping it an sp object and using "over" will be much faster than using "extract". You will also be processing everything in memory. – Jeffrey Evans Jan 27 '14 at 19:35
3

The example data from Jeffrey

library(raster)
r <- raster(ncols=10, nrows=10)
set.seed(0)
x <- runif(ncell(r))
x[round(runif(25,1,100),digits=0)] <- NA
r[] <- x
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1),  Polygons(list(Polygon(cds2)), 2)))
polys <- SpatialPolygonsDataFrame(polys, data.frame(ID=sapply(slot(polys, "polygons"), function(x) slot(x, "ID"))))

Now use extract

extract(r, polys, fun=function(x, ...) length(na.omit(x))/length(x))
#[1] 0.8333333 0.6666667

If you have many rasters, first use stack to combine them (if they have the same extent and resolution)

To get the actual polygon area you should not use the slot(i, 'area') approach. For planar data you can use rgeos::gArea(polys, byid=TRUE) For spherical data (lon/lat) you can use geosphere::areaPolygon

3

I am not sure if you want the ratio based on the "real value" of the polygon(s) areas or the areas of the cells intersecting them. Here is some example code that uses all cells intersecting the polygons (basically, ratio of NA cells to non-NA cells). It is a dummy example and you will need to write your own function.

    # Create some example data
    require(raster)
    require(sp)

    r <- raster(ncols=10, nrows=10)
      x <- runif(ncell(r))
        x[round(runif(25,1,100),digits=0)] <- NA
          r[] <- x
      cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
        cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
          polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                                   Polygons(list(Polygon(cds2)), 2)))
            polys <- SpatialPolygonsDataFrame(polys, data.frame(ID=sapply(slot(polys, "polygons"), 
                                              function(x) slot(x, "ID"))))
plot(r)
  plot(polys, add=TRUE)

You can use this code snippet to add an area column to your polygon data by extracting from the area slot. This could be used if you want to ratio using the "real" polygon area.

# Add area of polygon(s)
polys@data <- data.frame(polys@data, Area=sapply(slot(polys, 'polygons'), 
                         function(i) slot(i, 'area')))  

The most efficient, and considerable faster, alternative to for loops are "apply" like functions. There are a number of these available in R that are utilized for different object classes or data structures. In this case, since extract returns a list, we will use lapply (list apply). This is a way to apply a base or custom function to a list object. The the object class stored in the list is a vector, the function is quite straight forward. If you use extract on a brick or stack raster object the resulting objects stored in the list would be data.frame objects.

# On a single raster object, extract returns list object with stored vectors.                           
( vList <- extract(r, polys, na.rm=FALSE) )
  class(vList)

# Use lapply to apply function that calculates ratio of NA to non-NA values
#   wrapping lapply in unlist() collapses result into a vector  
aRatio <- function(x) { if(length(x[is.na(x)]) > 0) (length(x[is.na(x)]) / length(x[!is.na(x)])) else 0 }  
  ( vArea <- unlist( lapply(vList, FUN=aRatio ) ) )

# Assign ordered vector back to polygons
polys@data <- data.frame(polys@data, NAratio=vArea)
  str(polys@data)         
  • Thanks Jeffrey! I've learned quite a bit from your answer. But I think I didn't explain myself well enough. The ratio that I'm after is Area of NonNA Cells within Poly1 to Area of Poly1. Some polygons are not entirely covered by raster cells. Writing the mean value of all cells within a polygon into vList is great. Now I only need to get the number of NonNA cells as well from which the mean was derived as I know the area of each cell. The ratio can then be easily derived by (number of cells * cell area) / polygon area. Many Thanks! – Hubert Jan 28 '14 at 16:25
1

I do not have access to your files, but based on what you described, this should work:

library(raster)
mask_layer=shapefile(paste0(shapedir,"AOI.shp"))
original_raster=raster(paste0(template_raster_dir,"temp_raster_DecDeg250.tif"))
nonNA_raster=!is.na(original_raster)
masked_img=mask(nonNA_raster,mask_layer) #based on centroid location of cells
nonNA_count=cellStats(masked_img, sum)

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.