TransWikia.com

Map all time steps of a NetCDF file in R using a loop

Geographic Information Systems Asked by Maria Karypidou on September 9, 2020

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" ;

One Answer

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.

Answered by Sam on September 9, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP