I have a raster layer of approximately 240 million cells. It could be a DEM or anything similar. I also have a series of sampling-site polygons of irregular shape located throughout the extent of the DEM but only occupy about 0.1% of the space. I am using statistical tests to compare the distribution of values (elevations in this example) between the sample-site locations and the overall DEM (considered background values).

I have used the MASK() function in the Raster package of R to extract the the DEM values where they intersect non-Null values of the sample-site raster. It returns a raster that is 99% NA and the remainder are the values at the sample sites.

From here, I draw roughly 50K random samples from the DEM and I test them against the distribution of values at sampling-site locations. However, I have only found the most inefficient ways to remove or ignore the NA values (eg. x <- na.omit(y))

What is an efficient way to extract the values from the sampling-site raster into a vector that contains no NA values?

I aware of the na.rm flag fro some functions in the raster package, but I have not used it to success.
I have successfully implemented this general process by turning the sampling-site raster into cell centroids within sampling-sites and extracting values, but that has extra steps.

This is the most effective way I have found site_smpl_freq <- as.data.frame(freq(site_smpl, useNA='no', progress='text') but it still takes a few minutes to run

This is a code example.

    ###sample background data and put in data frame
    rand_smpl <- sampleRandom(raster("DEM"), collapse =''))),50000)
    rand_smpl_freq <- as.data.frame(count(rand_smpl))

    ### then a bit that modifies the table if values exceed a certain threshold

    ### extract raster values at the location of sampling-sites
    sites_extract <- mask(DEM, sampling_sites_raster)

    ### Here is where I create the freq table using the useNA option
    ### It works, but it takes a long time and seems inefficient
    site_smpl_freq <- as.data.frame(freq(sites_extract, useNA='no', progress='text'))

    ### then I modify the table as above if values exceed a certain threshold
    ### This is followed by some manipulations that add columns for cumulative counts and percentages, then to ggplot()

2 Answers 2


I really do not follow the logic of what you are attempting with "freq". The DEM values should be continuous and using freq just does not make sense to me but, you may have an unstated reason for doing so.

I would follow @Thomas advice and use extract. You will then have a list object with a vector of raster values for each sample site. You can use lapply to calculate moments to compare to a random sample of the entire raster.

Here is a quick example of one way to approach your problem.


# Create some example data
r <- raster(ncol=500, nrow=500)
  r[] <- runif(ncell(r),0,1)
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)))

  plot(polys, add=TRUE)                           

# Create a vector of raster values within polygons,
#   here I collapse the list into a vector
v <- unlist(extract(r, polys))

# Create random sample of raster
rs <- sampleRandom(r, length(v), na.rm=TRUE)

# Evaluate moments of each distribution

# Test distributional equlaity using Kolmogorov-Smirnov Test
ks.test(v, rs, alternative="two.sided", exact=FALSE)  

# Plot distributions for visual comparison 
rast.den <- density(rs)
samp.den <- density(v)
rast.den$y <- rast.den$y / max(c(rast.den$y, samp.den$y)) 
samp.den$y <- samp.den$y / max(c(rast.den$y, samp.den$y)) 
plot(rast.den, type="n", main="", 
  xlim=c(min(c(rast.den$x, samp.den$x)), max(c(rast.den$x, samp.den$x))), 
  ylim=c(min(c(rast.den$y, samp.den$y)), max(c(rast.den$y, samp.den$y))) )     
        polygon(rast.den, col=rgb(1,0,0,0.5))
        polygon(samp.den, col=rgb(0,0,1,0.5))
     legend("topright", legend=c("Raster","Sample"), 
            fill=c(rgb(1,0,0,0.5),rgb(0,0,1,0.5))  )

# FYI, here is how you apply a function on the results of extract,
#   accounting for NA's  
v <- extract(r, polys)
unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))

If you mask the raster by raster, you will always get another huge raster. I don't think this is a way to make things faster.

What I would do is to try to mask by polygon layer using extract:

res <- extract(raster, polygons)

Then you will have all the cell values for each polygon and can run freq on them.

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.