Hello, I am new to R so probably a basic question...
I am working with climate data that spans 1000 years in length in monthly increments. The data is a variable, say temperature, on a spatial grid with 144 longitudes and 90 latitudes (in other words, it spans the globe with 2.5 x 2 degree resolution). The dimensions are thus 144 x 90 x 12000, and for example, saying temp[10,19,345] would output the temperature at the 10th longitude grid box, 19th latitude, and 345th month.
Right now I am working on a program that produces a time-series of average temperatures over a prescribed region. That is, I want to define a "box" (say from the 10th to the 30th longitude, and from the 6th to the 8th latitude) and produce a vector with one column (of the average temperature) and 12,000 rows.
Right now, I have done the spatial averaging shown below (please assume I have already loaded in the data and declared names to the temperature and area variables used):
Also note that "area" variable is needed because all the grid boxes do not have the same surface area, and so I have to weight each box accordingly. The surface area for a given box changes with latitude but is independent of its longitude. The area of a given box is given as area[longitude,latitude].
######inputs######
# declare time range of interest
timecall1= 1
timecall2=12000
lat1=6 # latitude range to do average
lat2=7
lon1=5
lon2=6 # longitude range to do average
##################
k =lat1; sum=0;
while (k != lat2+1) {
sum=sum+sum(area[lon1,k]*temp[lon1:lon2,k,timecall1])
k= k+1;
}
avg_value=sum/(sum(area[lon1:lon2,lat1:lat2]))
Basically what I'm doing is averaging all the longitude values of temperature for a particular latitude, then moving up to the next latitude and repeating the process, etc.
I've checked the spatial averaging and it seems to be fine, though I come from a science and not a programming background, so if there is a more efficient way to do the above please let me know. I can imagine that with a data set this big such a calculation might get computationally expensive.
The last step is to make it a time series. The above program only works for a single month. It seems trivial but I can't get it. I've tried changing the value in the parentheses to "timecall1:timecall2" or adding in another loop, but it just ends up adding all the values in a way that I don't want.
My initial thought is to store this in R as a 3D array. It's like a matrix object, but it allows higher dimensions. In this case, you'll have spatial panels with another dimension of time. Each panel is a month. For prototyping, consider this
So what you're asking for is to take a subset of a panel and calculate its average for all 12,000 panels. The result can be stored simply as a 12,000 point vector in R, each position corresponding to its respective panel month. How would we simulate that in our prototype here?
Well, if I want to only look at rows 4 through 6, it's as easy as
Code:
a[4:6, , ]
I can of course include such subsets for the columns and the panel dimension. Both the syntax and concept here correspond to what you wanted. Of course, if we think of the rows in my array here as latitude and the columns as longitude, I don't have much to demonstrate! Still, suppose this was your data set and you wanted to look at the box you described. It would be as simple as
Code:
a[10:30, 6:8, ]
This would give us the box on each month. The question now is how do I take the average on each of the panels? Using my example above, we simply use apply on it with the MARGIN parameter on the panels
Code:
apply(a[4:6, , ], 3, mean)
This gives us
Code:
[1] 44.00000 73.83333 64.16667
If you took the 6 values on each of the panels of a[4:6, , ], you'll find that it matches those 3 vector point values. If you did this with the subset indexes for the 2 dimensions, you should end up with a vector 12,000 points long of the means you wanted.
The question you need to ask yourself now is how do you populate the array with the data you have. That, however, requires a bit more knowledge about how your data is stored and pulled into R. This is some pretty large data, but not the sort that require particularly special care. To save time, however, you will want to use R effectively. So, if you give us more details about your data files, we could prescribe a way to scan them into R to store them in a 3D array.
I can test this and provide further details tomorrow when I am at work.
For now, I can elaborate by saying I am using netCDF files and using "ncdf" in order to manipulate the data within the R environment.
For example, I say temp = get.var.ncdf(nc,"tsurf") for the "tsurf" variable contained within nc. nc=open.ncdf(file name.nc). I'm not sure if you are familiar with ncdf or are looking for something more fundamental, but I was not involved with the construction of the netCDF files. I'm only bringing them into R to do some data analysis, but I think you can treat them as any arbitrary 144 x 90 x 12000 matrix....
I used to work next to an atmospheric scientist, so I know what netCDF is, but I don't know how that library pulls it into R. If you can take a snippet of the data and print the output from dput on that data object, we can dget to load it up into R ourselves. You can also show the print out of str(...) that should detail the structure of the object.
Here is some output. Not sure if this is what you are looking for (note this is 1200 not 12,000 months in length, since I have 10 files split up into centuries)...if not, let me know what you need me to do
And I'd still like to follow up on the original method I had started off with, even if it's not suitable for a large dataset...just for my own education.
Okay. Don't ever do your original method. Done. You don't iterate over vectors in R. It's slow as hell. You use vectorized operations. Just compare the speed of these two commands
Code:
x <- rnorm(10000000)
sum(x) # Instant! Vectorized functiion
y <- 0
for(i in 1:length(x)) y <- y + x[i] # SLOOOOW! Don't use loops
That's an ugly data set, but the values are accessible. It looks to me like you'll want to first extract the values and form your matrix. I'd focus on one decade file at a time. I'll have to look at this data more deeply to figure out exactly how it links together. You might consider looking at some of the spatial packages (e.g., sp), as they may have methods for handling this sort of data better.
Reply With Quote
