I am not 100% sure if GIS is the correct stackexchange for this, but here we go.

A problem I have been bested by: 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))

*edited to add: 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. I apologize for being so verbose, but I would appreciate any and all ideas.

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

edit* to add 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()
  • A code example would be most welcome. – Jeffrey Evans Nov 22 '13 at 22:00

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 ))
  • Thank you both for addressing my query. Tomas has a good approach and Jeffrey gave a very good and reproducible example. Further, thanks for the plot example, that is very useful. My challenge now will be to a list of polygon coordinates for the shapefile of sampling locations. I will look to the shapefile package for some inspiration. PS - regarding the use of freq(), at times I can doing this routine to pull values from discreet rasters to compare the outcomes of various models. This will work for both continuous and discreet cases. – Mr.ecos Nov 27 '13 at 20:09
  • Update: I used the maptools package to get the sampling site shapefile into a spatial object. sampling_locals <- readShapePoly(site_Locals_SHP, IDvar=ID_column, proj4string=CRS(Prj)) then followed the logic you provided above with the exception of removing the NA values with var_sample <- var_sample[!is.na(var_sample)] instead of side-steeping them with na.rm=TRUE Thanks again! – Mr.ecos Dec 2 '13 at 17:09
  • I would recommend using readOGR in the rgdal package. This will retain the projection information. – Jeffrey Evans Dec 2 '13 at 17:15

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.