On NCDF Climate Datasets

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 <- 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),
+ month.abb[month(dates)] %>%
  tolower, sprintf("%02s", day(dates)),
  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.


13 thoughts on “On NCDF Climate Datasets”

  1. I wonder, will you be talking about grib and grib2 files as well? Those both seem to be more and more common, and there don’t appear to be as good resources for both.

  2. This didn’t work on my Windows 7 system. I couldn’t run the open.ncdf function, though the file is present and there were no obvious errors in the file download. Thoughts?

    > nc dir()

    [13] “temp” “temp.nc”
    > file.size(“temp.nc”)
    [1] 391738122

    1. I ran into this too. Even though the download command ‘works’, because you’re on a windows machine it needs an additional flag: mode=”wb”. In other words:

      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”,mode=”wb”, destfile=”temp.nc”)

      Hope this helps!

  3. Thank you for the tutorial – I haven’t worked with ncdf in a while so it is good to get a refresher on it. I tried to follow along the script, but can’t find the code to create the variables ind_conserver or dates_conserver. Perhaps I missed something, but, if necesssary, could you edit the post to add those lines?

Leave a Reply to Arthur Charpentier Cancel reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.