3

I have 20 rasters that cover the same geographic extent, are in the same projection, have the same cell size etc. Instead of a stack of rasters, I want a table where the columns are the different rasters and the rows are the pixels. I have tried multiple tools in ArcGIS, such as Extract Multi Values to Points and Sample, as well as using R. The issue is my rasters are very large, about 60 million pixels each. Those tools either don't work correctly and output incorrect values (Extract Multi Values to Points) or fail after a short time (Sample). I have tested them out with much smaller sub-samples (i.e. 100 points) and they work fine in that case. My R attempts have been met with errors due to lack of memory. Any tips or other ways I can approach this problem? I am willing to try different methods! I have read similar questions here and elsewhere in the past but it appears the large size of my dataset is leading to additional issues.

The goal is to have a table like this:

    • R1 - R2 - R3 etc.
  • P1
  • P2
  • P3 etc.

That I can bring into R and work with further.

EDIT: My R code... very simple at the moment (also down to 12 rasters).

Rasters <- stack(NDVI_C, Landcover, Elevation, Slope, Aspect, TPI, TWI, Lake_50, River_250, MajRiver, FTE, Precip, Temp)
d <- values(Rasters)
test <- values(NDVI_C)

A couple issues: memory and dealing with NA values. Want final table to only include pixel stacks that have no NA values (i.e. every Raster has a value for that pixel)

  • Library raster in R can read raster object without loading them into the memory. Because your rasters are big, you create a database on disk and fill it step by step while extracting raster values by blocks that can be handled in your memory with raster::getValuesBlock – Sébastien Rochette May 2 '17 at 17:24
  • @StatnMap okay thanks for the idea I will try to figure that out – Mitchell May 2 '17 at 19:26
  • To drop row NA's just use na.omit() or complete.cases() – Jeffrey Evans May 2 '17 at 22:28
3

This is, generally, a very bad idea and defeats the point of a raster format. There is very little in the way of software that can deal with a 60M x N dataset. With 20 parameters the number of observations grows to n=1,200,000,000. In statistical terms, you would be better served taking a sub-sampling approach or just leverage raster functions for your analysis.

To directly address your question it is easy to convert a raster stack to a data.frame object (in your terms, a table) that could be exported to whatever format you would like.

Here you can coerce a raster stack to a data.frame object. To write to a readable format for other software you can use write.csv().

library(raster)
r <- stack(system.file("external/rlogo.grd", package="raster")) 

r.df <- as.data.frame(getValues(r)) 
  head(r.df)

If you need a memory safe approach to create an ASCII file you can use read the raster stack in blocks and the use write.csv with the append = TRUE argument. If set up in a for loop you can literately stream data into an ASCII file.

b <- getValuesBlock(logo, row=35, nrows=3, col=50, ncols=3, lyrs=1:3)
write.csv(b, "raster.csv", col.names = TRUE, row.names = FALSE, append = TRUE)

A sampling approach can capture the sample variation and allow for statistical analysis without having to account for the entire population. Here is an example with the above data where we draw a n=1000 random sample and display the original and sample distributions.

r.df <- as.data.frame(getValues(r))
rs <- as.data.frame(sampleRandom(r, 1000))

par(mfrow=c(2,1))
  plot(density(r.df[,1]), main="Original distribution band 1")
  plot(density(rs[,1]), main="Sample distribution band 1")

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.