1

I have a NetCDF containing temperature values over a certain region for 2011, with 1-hour intervals, so that means that I have 8760 maps. What I want to do, is plot all 8760 maps using one color palette, and also on every map plot the exact day and time on the map, because I want to show the evolution of temperature over a region, in a year. How can I do this using a for loop in R?

The ncdump -h of my .nc files gives the following result:

netcdf temperature {
dimensions:
    rlon = 169 ;
    rlat = 155 ;
    height = 1 ;
    time = UNLIMITED ; // (8760 currently)
variables:
    double rlon(rlon) ;
            rlon:standard_name = "grid_longitude" ;
            rlon:long_name = "longitude in rotated pole grid" ;
            rlon:units = "degrees" ;
            rlon:axis = "X" ;
    double rlat(rlat) ;
            rlat:standard_name = "grid_latitude" ;
            rlat:long_name = "latitude in rotated pole grid" ;
            rlat:units = "degrees" ;
            rlat:axis = "Y" ;
    char rotated_pole ;
            rotated_pole:grid_mapping_name =     "rotated_latitude_longitude" ;
            rotated_pole:grid_north_pole_latitude = 39.25 ;
            rotated_pole:grid_north_pole_longitude = -162. ;
    double height(height) ;
            height:standard_name = "height" ;
            height:long_name = "height" ;
            height:units = "m" ;
            height:positive = "up" ;
            height:axis = "Z" ;
    double time(time) ;
            time:standard_name = "time" ;
            time:units = "hours since 2011-01-01 00:00:00" ;
            time:calendar = "proleptic_gregorian" ;
    float var11(time, height, rlat, rlon) ;
            var11:table = 2 ;
            var11:grid_mapping = "rotated_pole" ;
3
  • 3
    Start with library(raster); r <- raster("ncfile.nc", band = 1); plot(r); Use "zlim = "argument to plot to keep the range constant, and pass in the colours you want to "col = " – mdsumner Jan 22 '17 at 2:53
  • How could this work in a loop, printing each time the exact day and time of the map? – Maria Karypidou Jan 22 '17 at 13:42
  • Try getZ(r) does that give the right time? – mdsumner Jan 22 '17 at 21:39
0

This is given the assumption that there is only one variable in the netcdf, called temperature.

Check the name of your dimenions (assuming it's called "time");

names(your.nc$dim)

then use the length of your time dimension for the loop. It will loop through every time step, extract the raster and plot it. You'll need to determine what text you add etc, but the e.g. uses POSIXct notation:

require(ggplot2)
require(raster)
require(ncdf4)
# loop through length of timesteps (your case = 8760)

for(i in 1:length(your.nc$dim$time$vals)){

# rasterize data for every time step (for variable, "temperature")
  r <- raster("./path/your.nc",varname="temperature",band=i)
  rxy <- as.data.frame(r,xy=T)
# plot in ggplot2
  ggplot()+
    geom_tile(data=rxy,aes(x=x,y=y,fill=rxy[,3]))+
    labs(x="X",y="Y",title=as.POSIXct(i*3600,origin='2011-01-01 00:00'))

}

the as.POSXIct function takes a number, in seconds, and converts it to a datetime object from the origin nominated.

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.