R loop: name objectives by variable

Question

I need to read daily netcdf files every month and give each a name ended with the date

library(raster)
year<-2004
startmonth<-1
for(monthd in 31){
    days<-formatC(monthd, width=2, flag="0")
    month<-formatC(startmonth,width=2,flag="0")
    sm=raster(paste(year,month,days,"1.nc",sep=""),varname="sm")
    monthd<monthd+1
}

In the end I should have raster objectives named as sm01 sm02 . . . sm31

for January. There must be a simple way to do it, I'm just very fresh in coding.

Or better yet, use `list.files` and then stack the objects into a `rasterstack` using `stack` — Simon O'Hanlon, May 02 '13 at 20:24

score 1 · Accepted Answer · edited Jun 20 '20 at 09:12

You want to take a mean of a set of raster files. Package raster has built-in object types for handling situations such as this. It is much more efficient to create a rasterStack and use calc to find the mean:

days<-formatC(1:31, width=2, flag="0")
files <- list( paste("2004" , "01" , days , "1.nc" , sep="" ) )

## First 3 filenames:
files[[1]][1:3]
# [1] "200401011.nc" "200401021.nc" "200401031.nc"

sms <- stack( files )
smMean <- calc( sms , mean , na.rm = TRUE )

Edit based on OP comment

You cannot pass a list of filenames directly to raster. From the manual, it states:

x: filename (character), Extent, Raster*, SpatialPixels*, SpatialGrid*, object, 'image', matrix, im, or missing.

Therefore if you must have individual raster objects (and I strongly advise against it in this case) then you could use your loop as originally planned, or you could use:

smRaster <- lapply( files , raster , varname = "sm" )

This will return one list object, each element of which is a raster object. This is probably not that useful to you, but you can then access each one using smRaster[[1]] , smRaster[[2]] etc.

But use a raster stack if your files have the same extent and resolution!

Reading the files in using stack is likely to be more efficient and help you write more readable, shorter code. You can operate on all the rasters at once using convenient syntax, e.g. if I wanted to make a new raster that showed the sum of all the other rasters:

## First we use our files list to make the rasterStack
smStack <- stack( files , varname = "sm" )

## 'sum' returns a new rasterLayer, each cell is the sum of the corresponding cells in the stack
smSum <- sum( smStack , na.rm = TRUE )

## 'mean' does the same!!
smMean <- mean( smStack , na.rm = TRUE )

## What if we want the mean of only the first 7 days of rasters?
smMean <- mean( smStack[[1:7]] , na.rm = TRUE )

## To see how many rasters in the stack...
nlayers( smStack )

## To make a new object of just the first 7 rasters in the stack
subset( smStack , 1:7 )

## Is roughly equivalent to using the `[[` operator like this to get the first 7 layers in a new object
newRaster <- smStack[[1:7]]

## If you want to access an individual layer..
smStack[[ layernumber/or layername here ]]

## And to get the names:
names( smStack)

I would really strongly advise using a stack if you rasters cover the same spatial extent and resolution. It is more efficient, both for R and for your coding to just use a stack! I hope I have convinced you! :-)

I agree stack is a more elegant way doing this. However did you miss a step converting the netcdf files to raster? I found it difficult to do with `smRaster<-raster(files)` Any suggestions? — EDU, May 02 '13 at 23:51
@user2271576 I updated the question. In short you can't pass a list of files to `raster`, ***but*** don't use `raster` in this instance use `stack`!! I hope after reading the question I will have convinced you. I *think* it will really aid your work flow. HTH. — Simon O'Hanlon, May 03 '13 at 09:55
I've already been convinced by using stack, since I thought about stack in the first place, but wanted to avoid adding complexity to the process. I've tried to stack "files" by using `smStack <- stack( files , varname = "sm" )`. However it's giving 10 warnings with the first one "1: In if (x == "" | x == ".") { : the condition has length > 1 and only the first element will be used". Also the smStack only has one layer rather than the desired 31 layers. I've put my soil moisture data here: [link](https://www.dropbox.com/sh/47sor6faoi0bsfl/5o6w2g7JVa). — EDU, May 03 '13 at 17:10
@EDU I'll have a look into it tomorrow. Have you named your files correctly? If the files are all in the same directory, just use `files <- list.files( "/path/to/directory/" , pattern = ".nc" , full.names=TRUE)` — Simon O'Hanlon, May 03 '13 at 18:23
I think I've found the reason that stack couldn't read the list. We should use `list.files` rather than `list`, as `list.files` produces a character vector of the names whilst `list` produces a real 'list'. Obviously the `stack` recognizes character vector better. Thank you all very much. I think my question has been fully solve at this stage. Y'all have a pleasant weekend! — EDU, May 03 '13 at 18:55

score 0 · Answer 2 · answered May 02 '13 at 19:19

0

First thing: your for loop is off, you only iterate over one value, 31. You need to change it to:

for (monthd in 1 : 31) { …

Next, to create a variable based on a string, use assign:

assign(paste('sm', monthd, sep = ''), someValue)

Finally, remove the last line of your loop – it’s syntactically wrong and it achieves nothing, once corrected.

However, I strongly advise against using assign here. What you want isn’t lots of individual variables, you simply want a vector. You can simply create it like this:

sm <- vector(length = 31)
# or
sm <- 1 : 31
# or, same as previous
sm <- seq(1, 31)
# or
sm <- rep(someValue, 31)

answered May 02 '13 at 19:19

Konrad Rudolph

530,221
131
937
1,214

More robust : for ( monthd in length(n.days) ) where you've determined the number of days in the given month programmatically, e.g. [link] (http://stackoverflow.com/questions/6243088/find-out-the-number-of-days-of-a-month-in-r) – c.gutierrez May 02 '13 at 20:02
@c.gutierrez Now you’ve made the same error as OP, you need to replace `length` with `seq` or use `1 : n.days` – Konrad Rudolph May 02 '13 at 20:11
First, thanks for all the help. The reason I want to write 31 rasters is I'll later on do `sm_jan<-mean(sm1,sm2...sm31,na.rm=TRUE)`. How may I assign each raster to a list/vector such as sm with values 1 to 31? – EDU May 02 '13 at 20:54
@KonradRudolph my mistake; so becomes : for ( monthd in 1:n.days ) where n.days is returned value from function in above-mentioned post. In the below answer, line 1 becomes days<-formatC(1:n.days, width=2, flag="0" . – c.gutierrez May 02 '13 at 22:08
@user2271576 Not sure what you mean by the assignment, but your code becomes simply `sm_jan <- mean(sm, na.rm=TRUE)` and this works *in general*! – Konrad Rudolph May 02 '13 at 22:46
I think the gist of this should be to convince the OP to use a `stack` in this instance. It will be a *lot* more convenient for their use-case! – Simon O'Hanlon May 03 '13 at 09:57
@SimonO101 Maybe, I don’t know the `raster` package (don’t work with spatial data) and was just answering the general question about vector / loop usage. I trust your answer is correct though. – Konrad Rudolph May 03 '13 at 10:11

R loop: name objectives by variable

2 Answers2

Edit based on OP comment

But use a raster stack if your files have the same extent and resolution!