On NCDF Climate Datasets

September 3, 2015
By

(This article was first published on Freakonometrics » R-english, and kindly contributed to R-bloggers)

Mid november, a nice workshop on big data and environment will be organized, in Argentina,

We will talk a lot about climate models, and I wanted to play a little bit with those data, stored on http://dods.ipsl.jussieu.fr/mc2ipsl/.

Since Ewen (aka @3wen) has been working on those datasets recently, he kindly told me how to read those datasets (in some ncdf format). He did show me interesting functions, and did share with me some lines of codes. I will mention some of them in a brief post, because I’ve been impressed to see how memory issues were not a big deal, here… But first of all, we need half a dozen packages

> library(ncdf)
> library(lubridate)
> library(dplyr)
> library(stringr)
> library(tidyr)
> library(ggplot2)

Then i will download one of those files, from January 2000 till December 2012.

> download.file("http://dods.ipsl.jussieu.fr/
mc2ipsl/CMIP5/output1/IPSL/IPSL-CM5A-MR/historicalNat/day/atmos/day/r1i1p1/latest/tas/tas_day_IPSL-CM5A-MR_historicalNat_r1i1p1_20000101-20121231.nc",destfile="temp.nc")

Those file are not huge. But big.

> file.size("temp.nc")
[1] 390965452

That’s almost a 400Mo file. For one single file (and there are hundreds on those servers). Now, if we ‘read’ that file, in R, we have a rather small R object,

> nc <- open.ncdf("temp.nc")
> object.size(nc)
155344 bytes

That’s only 150Ko… One out of 2516 compared with the file online. Somehow, we only extract partial information here. We will scan the file completely later on… We can already create some variables that will be used later on, such as the date, the latitude, the longitude, etc

> lon <- get.var.ncdf(nc,"lon")
> (nlon <- dim(lon))
[1] 144
>  lat <- get.var.ncdf(nc,"lat")
> nlat <- dim(lat)
> time <- get.var.ncdf(nc,"time")
> tunits <- att.get.ncdf(nc, "time", "units")
> orig <- as.Date(ymd_hms(str_replace(
  tunits$value, "days since ", "")))
> dates <- orig+time
> ntime <- length(dates)

Here we have the range of latitude and longitude in the dataset, as well as the date. Now it is time to read the file. Based on the starting date, and asking only for temperature

> commencement <- ind_conserver[1]
> dname="tas"

we can red the file

> tab_val <- get.var.ncdf(nc, dname, start=c(1,1,commencement))
> dim(tab_val)
[1]  144  143 4745
> object.size(tab_val)
781672528 bytes

The file is big now, a bit less that 800Mo.

> fill_value <- att.get.ncdf(nc, dname, "_FillValue")
> tab_val[tab_val == fill_value$value] <- NA

Let us now summarize all the information in nice data frames, that will be used to generate graphs

> tab_val_vect <- as.vector(tab_val)
> tab_val_mat <- matrix(tab_val_vect, nrow = 
  nlon * nlat, ncol = ntime)
> lon_lat <- expand.grid(lon, lat)
> lon_lat <- lon_lat %>% mutate(Var1 = 
  as.vector(Var1), Var2 = as.vector(Var2))

Note that the file is still big (the same size as the previous file, here we just reshape the dataset)

> res <- data.frame(cbind(lon_lat, 
  tab_val_mat))
> object.size(res)
782496048 bytes

Here is the dimension of our dataset

> dim(res)
[1] 20592  4747

We can give understandable names to variables

> noms <- str_c(year(dates_conserver),
+ month.abb[month(dates_conserver)] %>%
  tolower, sprintf("%02s", day(dates_conserver)),
  sep = "_")  
> names(res) <- c("lon", "lat", noms)

Now, let us focus on the ploting part. If we focus on Europe, we need to filter the dataset

> res <- res %>%
  mutate(lon = ifelse(lon >= 180, yes = lon - 
  360, no = lon)) %>%
  filter(lon < 40, lon > -20, lat < 75, 
  lat > 30) %>% 
  gather(key, value, -lon, -lat) %>%
  separate(key, into = c("year", "month", 
  "day"), sep="_")

Let use filter, with only years from 1986 till 2005

> temp_moy_hist <- res %>%
+ filter(year >= 1986, year <= 2005) %>%
+ group_by(lon, lat, month) %>%
+ summarise(temp = mean(value, na.rm = TRUE)) %>%
+ ungroup

Now we can plot it. With the temperature in Europe, per month.

> ggplot(data = temp_moy_hist,
  aes(x = lon, y = lat, fill = temp)) + 
  geom_raster() + coord_quickmap() +
  facet_wrap(~month)

Actually, it is possible to add contour lines of European countries,

> download.file(
"http://freakonometrics.free.fr/carte_europe_df.rda",destfile="carte_europe_df.rda")
> load("carte_europe_df.rda")

If we aggregate all the month, and get only one graph, we have

> ggplot() +
  geom_raster(data = temp_moy_hist,
  aes(x = lon, y = lat, fill = temp)) +
  geom_polygon(data = carte_pays_europe, 
  aes(x=long, y=lat, group=group), fill = NA, 
  col = "black") +
  coord_quickmap(xlim = c(-20, 40), 
  ylim = c(30, 75))

I will spend some time now to analyze the maps, but I have been really impressed by those functions, that can partially read a ‘big’ dataset, to extract some information, and then we can read the file to create a dataset that can be used to get a graph.

To leave a comment for the author, please follow the link and comment on their blog: Freakonometrics » R-english.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)