WORKING WITH NetCDF IN R
MCS507 Lecture Note
GBODE Imoleayo Ezekiel
29/03/2025
1
#Introduction
• What is NetCDF?
– NetCDF – Network Common Data Form
– an array based data structure for storing multidimensional data
– written with an ASCII header and stores the data in a binary format.
##Advantages
• Saves space
– by writing binary files
– particularly because one does not need worry about the byte order.
– any byte- swapping is automatically handled by the netCDF libraries and;
• a netCDF binary file can thus be read on any platform.
• variety of data viewers already exist (e.g. ncview, R, NCL, CDO, FERRET, MATLAB, GraDs, IDL,
Data Explorer (OpenDX) etc.)
2
#Some features associated with the netCDF
• Coordinate systems: support for N-dimensional coordinate systems.
– X-coordinate (e.g., lat)
– Y-coordinate (e.g., lon)
– Z-coordinate (e.g., elevation/pressure)
– Time dimension
– Other dimensions
• Naming convention: classic NetCDF files are generally named with the .nc extension. + New version
of NetCDF is the NetCDF4 with .nc4 extension + High compressibility (i.e. about 20% less than .nc
file size)
• Variables: Support for multiple variables.
– E.g. precipitation, zonal (u) and meridional (v) wind components, temperature e.t.c.
• Geometry: Support for a variety of grid types (implicit or explicit).
– Regular grid (implicit)
– Irregular grid
– Points
• Self-Describing: Dataset can include information defining the data it contains.
– Units (e.g., mm, ms-1 , K or o C,. . . )
– Comments (e.g., titles, conventions used(e.g. CF, COARDS etc.), names of variables (e.g., zonal
wind, precipitation), names of coordinates (e.g., rectilinear, curvilinear, Gaussian coords..)
– Global attributes
See below links for further documentations on NetCDF data format
1. unidata NetCDF
2. unidata NetCDF Doc
3
File Structure
• Global Attributes: Describe the contents of the file.
• Dimensions: Define the structure of the data (e.g., time, depth, lat, lon).
– A netCDF dimension has both a name and a length.
– When a variable is defined, its shape is specified as a list of dimensions
– The number of dimensions is called the rank (dimensionality).
∗ a scalar variable has rank 0,
∗ a vector has rank 1 and;
∗ a matrix has rank 2.
• Variables: Holds the data in arrays shaped by dimensions
• Variable Attributes: Describes the content of each variable.
4
#Clear workspace
rm(list = ls())
########################################
#Install the necessary libraries
#1. [Link](c("PCICt", "ncdf4",
# "[Link]", "lubridate", "zoo",
# "fields","RColorBrewer", "maptools"))
#2. [Link]("maptools", repos = "[Link]
# 3. tinytex::install_tinytex()
#####################################################
#[Link]("chron")
#[Link]("sos")
########################################
library(ncdf4)
library([Link])
########################################
#Set working directory
#setwd("/Users/gbode/Documents/mcs/MCS507")
setwd("/Users/gbode/Documents/mcs/MCS507/class2024_2025")
######## load data
nc3 <- nc_open("/Users/gbode/Documents/mcs/MCS507/ERA5SfcTMinMax_mon_19810101-[Link]")
nc3 #ncdump data
## File /Users/gbode/Documents/mcs/MCS507/ERA5SfcTMinMax_mon_19810101-[Link] (NC_FORMAT_64BIT):
##
## 5 variables (excluding dimension variables):
## double time_bnds[bnds,time]
## float lon[longitude]
## standard_name: longitude
5
## long_name: longitude
## units: degrees_east
## axis: X
## float lat[latitude]
## standard_name: latitude
## long_name: latitude
## units: degrees_north
## axis: Y
## short tmax[longitude,latitude,time]
## long_name: Maximum temperature at 2 metres since previous post-processing
## units: K
## add_offset: 283.981992320285
## scale_factor: 0.00121006157896177
## _FillValue: -32767
## missing_value: -32767
## cell_methods: time: mean
## short tmin[longitude,latitude,time]
## long_name: Minimum temperature at 2 metres since previous post-processing
## units: K
## add_offset: 283.573540565479
## scale_factor: 0.00121903993957779
## _FillValue: -32767
## missing_value: -32767
## cell_methods: time: mean
##
## 4 dimensions:
## time Size:468 *** is unlimited ***
## standard_name: time
## long_name: time
## bounds: time_bnds
## units: hours since 1900-01-01 00:00:00.0
## calendar: gregorian
## axis: T
6
## bnds Size:2 (no dimvar)
## longitude Size:53 (no dimvar)
## latitude Size:53 (no dimvar)
##
## 6 global attributes:
## CDI: Climate Data Interface version 1.9.2 ([Link]
## Conventions: CF-1.6
## history: Mon Mar 1 11:05:45 2021: ncrename -v longitude,lon -v latitude,lat ERA5SfcTMinMax_m
## Mon Mar 01 10:58:17 2021: cdo monmean -sellonlatbox,2,15,2,15 /Users/gbode/Documents/climateAnalysis/
## Sat Aug 08 01:43:51 2020: cdo seldate,1981-01-01,2019-12-31 ERA5SfcTMinMax_day_19790101-[Link] E
## Sat Aug 8 00:44:41 2020: ncrename -v mn2t,tmin -v mx2t,tmax ERA5SfcTMinMax_day_19790101-[Link]
## Sat Aug 8 00:41:35 2020: ncrcat /Users/gbode/tutorial4python/era5/era5forAfrica/ERA5SfcTMinMax_day_1
## Mon Jul 27 10:06:49 2020: cdo daymean ERA5SfcTMinMax_1979.nc ERA5SfcTMinMax_day_1979.nc
## 2020-07-27 07:52:32 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /
## frequency: mon
## NCO: netCDF Operators version 4.9.5 (Homepage = [Link] Code = [Link]
## CDO: Climate Data Operators version 1.9.2 ([Link]
7
###read data variables
tmax <- ncvar_get( nc3, "tmax" )
lon <- ncvar_get( nc3, "lon" )
lat <- ncvar_get( nc3, "lat" )
tim <- ncvar_get( nc3, "time" )
tims <- [Link](nc3) #This function requires "[Link]"
#package
head(tims)
## [1] "1981-01-16 11:00:00" "1981-02-14 23:00:00" "1981-03-16 11:00:00"
## [4] "1981-04-15 23:00:00" "1981-05-16 11:00:00" "1981-06-15 23:00:00"
nc_close(nc3) ##close data file
##print variable dimension
dim(tim); dim(lon); dim(lat); dim(tmax) #note the data structure
## [1] 468
## [1] 53
## [1] 53
## [1] 53 53 468
#tmax(lon, lat, time)
##print lon and lat
"LON"; print(lon); "LAT"; print(lat)
## [1] "LON"
## [1] 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75
## [13] 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75 7.00 7.25 7.50 7.75
8
## [25] 8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.00 10.25 10.50 10.75
## [37] 11.00 11.25 11.50 11.75 12.00 12.25 12.50 12.75 13.00 13.25 13.50 13.75
## [49] 14.00 14.25 14.50 14.75 15.00
## [1] "LAT"
## [1] 15.00 14.75 14.50 14.25 14.00 13.75 13.50 13.25 13.00 12.75 12.50 12.25
## [13] 12.00 11.75 11.50 11.25 11.00 10.75 10.50 10.25 10.00 9.75 9.50 9.25
## [25] 9.00 8.75 8.50 8.25 8.00 7.75 7.50 7.25 7.00 6.75 6.50 6.25
## [37] 6.00 5.75 5.50 5.25 5.00 4.75 4.50 4.25 4.00 3.75 3.50 3.25
## [49] 3.00 2.75 2.50 2.25 2.00
9
##restructure the data
lon2 = lon#-180;
lat2 = rev(lat)
##print new LON and LAT
"LON2"; print(lon2); "LAT2"; print(lat2)
## [1] "LON2"
## [1] 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75
## [13] 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75 7.00 7.25 7.50 7.75
## [25] 8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.00 10.25 10.50 10.75
## [37] 11.00 11.25 11.50 11.75 12.00 12.25 12.50 12.75 13.00 13.25 13.50 13.75
## [49] 14.00 14.25 14.50 14.75 15.00
## [1] "LAT2"
## [1] 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75
## [13] 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75 7.00 7.25 7.50 7.75
## [25] 8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.00 10.25 10.50 10.75
## [37] 11.00 11.25 11.50 11.75 12.00 12.25 12.50 12.75 13.00 13.25 13.50 13.75
## [49] 14.00 14.25 14.50 14.75 15.00
tmaxN <- tmax[ , ncol(tmax):1 , ] #Reverse the lat/col
##Define longitude and latitude corners for Nigeria
maxlon <- 15
minlon <- 2
maxlat <- 15
minlat <- 2
##To subset time
library(PCICt)
library(lubridate)
10
##
## Attaching package: ’lubridate’
## The following objects are masked from ’package:base’:
##
## date, intersect, setdiff, union
mintime <- [Link]("1981-01-01", cal="proleptic_gregorian")
maxtime <- [Link]("2019-12-31", cal="proleptic_gregorian")
##Subset lon, lat and time for a region
idxlon <- which( lon2>= minlon & lon2 <= maxlon )
idxlat <- which( lat2>= minlat & lat2 <= maxlat )
idxtime <- which( tims>= mintime &
tims <= maxtime)
##Subset lon, lat and time for a specific site
slon = 7; slat = 7
idxlonS <- which( abs(lon2 - slon) == min(abs(lon2 - slon)))
idxlatS <- which( abs(lat2 - slat) == min(abs(lat2 - slat)))
#--split time
year = format(tims, format = "%Y")
month = format(tims, format = "%m")
day = format(tims, format = "%d")
hour = format(tims, format = "%H")
minute = format(tims, format = "%M")
second = format(tims, format = "%S")
idxtimeM <- which( tims>= mintime &
tims <= maxtime & month(tims)==5) # Extracting for only May
11
## Warning: tz(): Don’t know how to compute timezone for object of class PCICt;
## returning "UTC".
#--subset for a region
Tmax <- tmaxN[idxlon, idxlat, ]
llon <- lon2[idxlon]
llat <- lat2[idxlat]
tims2 = tims[idxtime]
tims3 = tims[idxtimeM]
#--subset for a specific site
TmaxS <- tmaxN[idxlonS, idxlatS, idxtime]
tmax.n <- tmaxN[idxlon, idxlat, idxtime]
dim(tmax.n)
## [1] 53 53 468
[Link] <- apply(tmax.n, 3, mean) #--domain average
tmax.mn2 <- apply(tmax.n, c(1,2), mean) #--time average
dim(tmax.mn2)
## [1] 53 53
#--aggregate data by time
[Link] <- aggregate([Link], by=list(year(tims2)), FUN="mean") #--year mean
## Warning: tz(): Don’t know how to compute timezone for object of class PCICt;
## returning "UTC".
[Link] <- aggregate([Link], by=list(month(tims2)), FUN="mean") #--monthly mean
## Warning: tz(): Don’t know how to compute timezone for object of class PCICt;
## returning "UTC".
12
"Nigeria LON";print(llon);
## [1] "Nigeria LON"
## [1] 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75
## [13] 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75 7.00 7.25 7.50 7.75
## [25] 8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.00 10.25 10.50 10.75
## [37] 11.00 11.25 11.50 11.75 12.00 12.25 12.50 12.75 13.00 13.25 13.50 13.75
## [49] 14.00 14.25 14.50 14.75 15.00
"Nigeria LAT"; print(llat)
## [1] "Nigeria LAT"
## [1] 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75
## [13] 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75 7.00 7.25 7.50 7.75
## [25] 8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.00 10.25 10.50 10.75
## [37] 11.00 11.25 11.50 11.75 12.00 12.25 12.50 12.75 13.00 13.25 13.50 13.75
## [49] 14.00 14.25 14.50 14.75 15.00
##Annual time series plot
plot([Link])
13
306
304
302
x
300
298
1980 1990 2000 2010 2020
Group.1
plot([Link], type = 'l')
14
306
304
302
x
300
298
1980 1990 2000 2010 2020
Group.1
plot([Link], type = 'l', xlab='Year', ylab='Maximum Temperature (degC)')
15
306
Maximum Temperature (degC)
304
302
300
298
1980 1990 2000 2010 2020
Year
#format ylab
plot([Link], type = 'l', xlab='Year', ylab=expression('Maximum Temperature ('*~degree*C*')'))
16
306
Maximum Temperature ( °C)
304
302
300
298
1980 1990 2000 2010 2020
Year
#add trend line to the time series plot
plot([Link], type = 'l', col = 'red', xlab='Year',
ylab=expression('Maximum Temperature ('*~degree*C*')'))
lm.p = lm([Link]$x ~ [Link]$Group.1)
abline(lm.p, col = 'grey', lwd = 3, lty=2)
17
306
Maximum Temperature ( °C)
304
302
300
298
1980 1990 2000 2010 2020
Year
##Annual/Seasonal cycle
plot([Link], type = 'l')
18
302
301
x
300
299
2 4 6 8 10 12
Group.1
plot([Link], col='red', type = 'l', xlab='Month',
ylab=expression('Maximum Temperature ('*~degree*C*')'))
19
302
Maximum Temperature ( °C)
301
300
299
2 4 6 8 10 12
Month
#customize x axes labels
plot([Link], type = 'l', xlab='Month', col='red',
ylab=expression('Maximum Temperature ('*~degree*C*')'),
axes=FALSE)
par(las=2)
labels = c('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul',
'Aug', 'Sep', 'Oct', 'Nov', 'Dec')
axis(1, at=1:12, labels = labels)
axis(2)
box()
20
Maximum Temperature ( °C)
302
301
300
299
Jan
Feb
Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
Nov
Dec
Month
library(fields) #for [Link] function
## Loading required package: spam
## Spam version 2.11-1 (2025-01-20) is loaded.
## Type ’help( Spam)’ or ’demo( spam)’ for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix ’.spam’ to the function name, e.g. ’help( [Link])’.
##
## Attaching package: ’spam’
## The following objects are masked from ’package:base’:
##
## backsolve, forwardsolve
21
## Loading required package: viridisLite
##
## Try help(fields) to get started.
library(RColorBrewer)
range(Tmax)
## [1] 286.4130 313.0029
[Link](llon, llat, Tmax[ , ,1], col = (rev([Link](10,"RdBu"))),
zlim=c(range(Tmax)))
14
310
12
305
10
llat
300
8
295
6
290
4
285
2
2 4 6 8 10 12 14
llon
#Compute the time average
TmaxN <- apply(Tmax, c(1,2), mean)
[Link](llon, llat, TmaxN)
22
14
302
12
300
10
298
llat
296
6
294
4
2
2 4 6 8 10 12 14
llon
#Compute the time average
TmaxN <- apply(Tmax, c(1,2), mean)
[Link](llon, llat, TmaxN)
#Add basemap
library(maptools)
## Loading required package: sp
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, were retired in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## [Link]
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## Please note that ’maptools’ will be retired during October 2023,
23
## plan transition at your earliest convenience (see
## [Link] and earlier blogs
## for guidance);some functionality will be moved to ’sp’.
## Checking rgeos availability: FALSE
data(wrld_simpl)
plot(wrld_simpl,add=TRUE)
14
302
12
300
10
298
llat
296
6
294
4
2
2 4 6 8 10 12 14
llon
#Make the plot appear better
image(llon, llat, TmaxN, col=rev([Link](1000, alpha = 1)),
xlab="Longitudes", ylab="Latitudes", main="Temperature Celsius")
contour(llon, llat, TmaxN, add=T, col= "grey", nlevel=20)
plot(wrld_simpl,add=TRUE)
24
Temperature Celsius
301 301
302.5
14
302
302.5
.5
300 301
300.5 301.
12
301 5
302 298.5 299
0.5 9
29
298
10
300.5 297 30
Latitudes
299.5 302
302
299 300
8
301 301.5
296 297 296
295
300 298
.5 3.5
6
8 29
29 297.5 8
29
7.5
296.5
29
299.5
297.5
4
297 297
9
29 6.5
29
297.5
2
2 4 6 8 10 12 14
Longitudes
#To save the plot in a device
png("ERA5_Tmax_Nigeria.png") #Open png device; can be pdf, jpeg, postscript etc
#Compute the time average
TmaxN <- apply(Tmax, c(1,2), mean)
[Link](llon, llat, TmaxN)
#Add basemap
data(wrld_simpl)
plot(wrld_simpl,add=TRUE)
[Link]() #Close the devise
## pdf
## 2
25
###########Panel Plot
# Here is a more learned strategy to add a common legend to a panel of
# plots consult the [Link] help file for more explanations.
# For this example we draw two
# images top and bottom and add a single legend color bar on the right side
# first divide screen into the figure region (left) and legend region (right)
[Link]( rbind(c(0, .8,0,1), c(.8,1,0,1)))
## [1] 1 2
# now subdivide up the figure region into two parts
[Link](c(2,1), screen=1)-> ind
# first image
screen(ind[1])
par(mar = c(1, 2, 0, 0))
image(lon2, lat2, tmaxN[ , ,1], col = (rev([Link](10,"Spectral"))),
zlim=c(range(tmaxN)), xlab="", ylab = "", xaxt = "n",
main = "Figure 1", useRaster = TRUE)
contour(lon2, lat2, tmaxN[ , ,1], col = 'red', add=TRUE)
plot(wrld_simpl,add=TRUE)
# Second image
screen( ind[2])
par(mar = c(2, 2, 0, 0))
image(lon2, lat2, tmaxN[ , ,1], col = (rev([Link](10,"Spectral"))),
zlim=c(range(tmaxN)),xlab="", ylab = "",
main = "Figure 2")
plot(wrld_simpl,add=TRUE)
?[Link]
# move to skinny region on right and draw the legend strip
26
screen(2)
[Link]( zlim=c(range(tmaxN)), [Link]=TRUE, smallplot=c(.1,.2, .3,.7),
col=(rev([Link](10,"Spectral"))))
Figure 1288
14
292 290
294
10
296
296
298 296
8
300
294
296
2
6
29
296 310
4
294
305
2
Figure 2 300
14
295
290
8 10
285
6
4
2
2 4 6 8 10 12 14
[Link]( all=TRUE)
27
#Practice Session
• Assignment 1: Plot blackbody emission (B_lambda) as a function of wavelength (lambda) for T=
300K and T= 6000K. Use the lambda range 0.1 - 100 micrometer). [To be submitted on Monday, 8th
March 2021]
• Assignment 2: Download the monthly sum rainfall data from CHIRPS ‘v2p0chirps_25_monSum.nc’
and reproduce similar plot in the ‘Working with NetCDF in R’ file fro the period 1986 - 2015 [To be
submitted on Thursday, 11th March 2021]
28