Read netCDF file with in R

Understand the metadata

We first explore the metadata of the file with the with the NetCDF() function from ncdump package.

metadata = ncdump::NetCDF("E:/MatlabWorking/Altimetry/old staff/wio_geostrophic_uv_july_2015.nc")

In summary, the metadata contains information about the netCDF file. One group of the data in the file include vector format variables like longitude, latitude and time as vector (table 1)

The second group of the data in the metadata are array of velocity. The array contained the zonal velocity (U) and meridional velocity (V) as matrix (table 2). Because this is the daily data, each matrice represent a single day.We can view the information about the file by printing the uv file in the console.

Extract the variable

Once we have identified the variables contained in the netCDF file, we use the nc_open() function to read the the file and assign it with uv name.

uv = nc_open("E:/MatlabWorking/Altimetry/old staff/wio_geostrophic_uv_july_2015.nc")
ssh = nc_open("E:/MatlabWorking/Altimetry/old staff/wio_ssh_july_2015.nc")

Once we have the file in the console, we are ready to extract the vector and array from the file. We can do this with ncvar_get()function.

## spatial components
lon = ncvar_get(uv, "lon")
lat = ncvar_get(uv, "lat")
## temporal component
time = ncvar_get(uv, "time")

## geogstrophic current
u = ncvar_get(uv, "u")
v = ncvar_get(uv, "v")

## sea surface height
adt = ncvar_get(ssh, "adt")

Printing the time we realized that time is the julian days but we do know the starting of the date. However, in the metadata provide information of the beginning data that we can use to transform this julian day into gregorian calendar that we are familiar with.

time

 [1] 23922 23923 23924 23925 23926 23927 23928 23929 23930 23931 23932
[12] 23933 23934 23935 23936 23937 23938 23939 23940 23941 23942 23943
[23] 23944 23945 23946 23947 23948 23949 23950 23951 23952

We can transform this days into the calender once we know the original date that AVISO used. Looking on the metadata we spotted that the original date for the calender was assigned as julian_day_unit: days since 1950-01-01 00:00:00. Here comes another challenges, the time is in Julian but the original date is in gregorian format. Therefore, we need to standardize the time to a common format. We first convert the origin time—1950-01-01 00:00:00 to julian day. Once the original time is in the same with the vector time file—julian format, we can add them up and convert from the julian to gregorian calender.

# convert time original (to) to julian 
to = insol::JDymd(year = 1950, month = 1, day = 1)

# add the original time to the extracted time
jd = to+time

#convert the julian day to gregorian calender
date = insol::JD(jd, inverse = TRUE)

The converted julian days are summarized in table 3 which show the julian number, julian day and the calender (gregorian date). We notice that our data were acquired from 2015-07-01 15:00:00 to 2015-07-31 15:00:00.

Transform to data frame

The final step is the transformation of matrix in arrays into data frame. Because there are thirty one matrix in each array—u,v, and adt, we chained the process in the loop as shown in the chunk below. I will not repeate the looping process in this post, If you cant follow the code in the chunk, I recommend you to read previous post that have sufficient information about how to loop and iterate repetitive process.

## zonal component

u.df = NULL

for (i in 1:length(date)) {

 ua = data.frame(lon,u[,,i] %>% as.data.frame())%>% 
   as.tibble() %>% 
   gather(key = "key", value = "u", 2:421) %>% 
   mutate(time = date[i], lat = rep(lat, each = 401)) %>% select(time, lon,lat, u)
 
 u.df = u.df %>% bind_rows(ua)

}


## meridional compoent

v.df = NULL

for (i in 1:length(date)) {

 va = data.frame(lon, v[,,i] %>% as.data.frame())%>% 
   as.tibble() %>% 
   gather(key = "key", value = "v", 2:421) %>% 
   mutate(time = date[i], lat = rep(lat, each = 401)) %>% select(time, lon,lat, v)
 
 v.df = v.df %>% bind_rows(va)
}


## sea surface height
adt.df = NULL

for (i in 1:length(date)) {

 adta = data.frame(lon,adt[,,i] %>% as.data.frame())%>% 
   as.tibble() %>% 
   gather(key = "key", value = "adt", 2:421) %>% 
   mutate(time = date[i], lat = rep(lat, each = 401)) %>% select(time, lon,lat, adt)
 
 adt.df = adt.df %>% bind_rows(adta)
}

Stitching the data

Once the geostrophic current (zonala and meridional) and the sea surface height anomaly data frame have been created, was combined and organized in a consistency format that makes analysis and plotting easy.

aviso = data.frame(u.df, v.df, adt.df) %>% 
  select(time, lon,lat, u, v, adt) %>% 
  as.tibble() %>% 
  mutate(day = lubridate::yday(time))

aviso$time = as.Date(aviso$time)
aviso$day = as.integer(aviso$day)

Table 4 is the random sample of twelve observations showing the information of geostrophic current and sea surface height within the tropical Indian Ocean for thirty one days in July 2015.

Visualize—Mapping the SSH

We have come long way from reading the data from the local machine into R’s environment, extract the variables and transform them into data frame. But the subtle information we are looking from this dataset are still hidden and we ought to uncover them. One way of understanding the data is through visualization. Visualizing data may range from common static exploratory data analysis plots to dynamic, interactive data visualizations in web browsers. R offers control over many aesthetic aspects of plots, but we will stick on ggplot2 developed by Hadley Wickham (2016) because it provide new ways to visualize and communicate data. Figure 1 for example show the sea surface height anomaly in the tropical Indian Ocean. We can easily spot region with higher and lower sea surface height anomaly. We see a region with high sea surface anomaly in the Mozambique channel and low sea surface height anomaly between latitude 8\(^\circ\)S and 4\(^\circ\)S and longitude 488\(^\circ\)E and 548\(^\circ\)E.

ggplot()+
  geom_raster(data = aviso %>% 
                filter(day == 190 & between(lon,38.5,54.5) & between(lat,-15,-.5)), 
              aes(x = lon, y = lat, fill = adt), interpolate = TRUE)+
  geom_sf(data = spData::world, fill = "grey85", col = 1)+
  coord_sf(xlim = c(39,53.5), ylim = c(-14.5,-1.5))+
  scale_fill_gradientn(name = "SSH\n(m)", colours = oce::oceColorsPalette(120))+
  labs(title = "", x = "", y = "")+
  theme_bw()+
  theme(axis.text = element_text(size = 11, colour = 1),
        legend.key.height = unit(1.5, "cm"))

Figure 1: Sea Surface Height Anomaly in the Tropical Indian Ocean Region as of 2015-07-01

If we want to check how the sea surface height anaomaly changes over the month of July 2015, then we ought to produce figure 1 for each day and obtain a total of 31 maps. Then visualize one after the other. That also wont save much and we will still miss the information as they are hidden in the table. Beside, we can animate them and make smooth transition of these maps. Thanks to Thomas Pedersen and David Robinson (2017)for developing a gganimate that transorm static plot of ggplot2 into animation. With few line of codes, we transformed static map in figure 1 into animated figure 2, which shown a seamless change of sea surface height in the area.

require(gganimate)

fig2 = ggplot()+
  geom_raster(data = aviso %>% filter(between(lon,38.5,60.5) & between(lat,-23,2)), 
              aes(x = lon, y = lat, fill = adt), interpolate = TRUE)+
  geom_sf(data = spData::world, fill = "grey85", col = 1)+
  coord_sf(xlim = c(39.5,58), ylim = c(-21.5,0.5))+
  scale_fill_gradientn(name = "SSH\n(m)", colours = oce::oceColorsPalette(120))+
  labs(title = "", x = "", y = "")+
  theme_bw()+
  theme(axis.text = element_text(size = 11, colour = 1),
        legend.key.height = unit(1.5, "cm"))+
  labs(title = 'Date: {frame_time}') +
  transition_time(time) +
  ease_aes('linear')
  
 animate(fig2)

Figure 2: Animation of Sea Surface Height Anomaly in the tropical Indian Ocean from July 1 to July 31, 2015

There are times we need to overlay geostrophic field that show speed and direction on the sea surfae height anomaly. Unfortunately, the density of the geostrophic often make the plot messy. Therefore, in the next post, we will reduce the density by creating equal size grid and then calculate the average velocity of the zonal (U) and meridional (V) velocity in each grid.