Access Oceanographic and Environmental Data with R

Geostrophic Current

One main advantage of using ERDDAP is that you only need to download the subset of the data you desire, rather than the entire dataset. You can subset based on temporal— time range and spatial—geographical extent of interest. This approach is convinient and essential for oceanographic and environmental variables, which are usually large.

Aviso Zonal Geostrophic Current is inferred from Sea Surface Height Deviation, climatological dynamic height, and ba. I will download the AVISO geostrophic current from erdTAgeo1day_LonPM180 servers, which contain two gridded Variables:

• u_current (Eastward Sea Water Velocity, ms^-1)
• v_current (Northward Sea Water Velocity, ms^-1)

Since we know the datasetID is erdTAgeo1day_LonPM180 and unsure whether is gridded or tabular dataset, you can check the metadata of the dataset with rerddap function info().

## obtain metadata of the dataset
geoinfo = info("erdTAgeo1day_LonPM180")
geoinfo

<ERDDAP info> erdTAgeo1day_LonPM180 
 Dimensions (range):  
     time: (1992-10-14T12:00:00Z, 2012-12-08T12:00:00Z) 
     altitude: (0.0, 0.0) 
     latitude: (-90.0, 90.0) 
     longitude: (-179.7625173611111, 179.9875) 
 Variables:  
     u_current: 
         Units: m s-1 
     v_current: 
         Units: m s-1 

If you do know the datasetId, the first way to find a dataset is to browse the builtin web page for a particular ERDDAP server. You can obtain a list of some public available ERDDAP servers with the servers() function as in the command below:

## call available ERDDAP servers
servers()

                                                                                            name
1                                                          Marine Domain Awareness (MDA) - Italy
2                                                                     Marine Institute - Ireland
3                                                       CoastWatch Caribbean/Gulf of Mexico Node
4                                                                     CoastWatch West Coast Node
5                     NOAA IOOS CeNCOOS (Central and Northern California Ocean Observing System)
6  NOAA IOOS NERACOOS (Northeastern Regional Association of Coastal and Ocean Observing Systems)
7                                         NOAA IOOS NGDAC (National Glider Data Assembly Center)
8    NOAA IOOS PacIOOS (Pacific Islands Ocean Observing System) at the University of Hawaii (UH)
9                     NOAA IOOS SECOORA (Southeast Coastal Ocean Observing Regional Association)
10                            NOAA NCEI (National Centers for Environmental Information) / NCDDC
11                                                NOAA OSMC (Observing System Monitoring Center)
12                                                           NOAA UAF (Unified Access Framework)
13                                                                   ONC (Ocean Networks Canada)
14                    UC Davis BML (University of California at Davis, Bodega Marine Laboratory)
15                                                                            R.Tech Engineering
16                                     French Research Institute for the Exploitation of the Sea
                                        url
1  https://bluehub.jrc.ec.europa.eu/erddap/
2           http://erddap.marine.ie/erddap/
3       http://cwcgom.aoml.noaa.gov/erddap/
4  https://coastwatch.pfeg.noaa.gov/erddap/
5     http://erddap.axiomalaska.com/erddap/
6           http://www.neracoos.org/erddap/
7       http://data.ioos.us/gliders/erddap/
8       http://oos.soest.hawaii.edu/erddap/
9            http://129.252.139.124/erddap/
10   http://ecowatch.ncddc.noaa.gov/erddap/
11             http://osmc.noaa.gov/erddap/
12     https://upwell.pfeg.noaa.gov/erddap/
13           http://dap.onc.uvic.ca/erddap/
14    http://bmlsc.ucdavis.edu:8080/erddap/
15            http://meteo.rtech.fr/erddap/
16  http://www.ifremer.fr/erddap/index.html

Surface Geostrophic Currents From Altimetry

The geostrophic equations are widely used in oceanography to estimate currents. The basic idea is to use hydrographic measurements of temperature, salinity or conductivitiy, and pressure to calculate the density filed of the ocean using the equation of state of sea water. There are several satellite-altimeter systems measuring the oceanic geostrophic current. I will not cover them in this post.

Briefly Theory of Geostrophic

Satellite altimeters used to measure surface geostrophic currents also measure wave height. The geostrophic balance requires that the Coriolis force balance the horizontal pressure gradient. The equations for geographic balance are derived from the equations of motion assuming the flow has no accelaration, \(du/dt = dv/dt = dw/dt = 0\): that horizontal velocities are much larger than vertical, \(w < u,v\); that the only external force is gravity; and that friction is small. With these assumptions we obtain equation (1)

\[ \begin{align*} \frac{\delta p}{\delta x} = \rho f v \\ \frac{\delta p}{\delta x} = \rho f u \\ \frac{\delta p}{\delta x} = -\rho g \:\tag{1} \end{align*} \]

Where \(2\Omega sin \varphi\) is the Coriolis pparamer. The equation can be writen as equation (2)

\[ \begin{align*} u = \frac{1}{f \rho} \frac{\delta p}{\delta y}\\ v = \frac{1}{f \rho} \frac{\delta p}{\delta x} \tag{2} \end{align*} \]

The geostrophic approximation applied at \(z = 0\) leads to a very simple relation— surface geostrophic current are proportional to surface slope. For example, consider a level surface slightly below the se sea surface, say two meters below the sea surface , at \(z = -r\) equation (3)

\[ \begin{equation} u = v = w = 0 \tag{3} \end{equation} \]

The pressure of the level is obtained with equation (4)

\[ \begin{equation} p = \rho g\; (\zeta + r) \tag{4} \end{equation} \]

Assuming \(\rho\) and \(g\) are essentially constant in the upper few meters of the ocean. Sustituting this into equation (1) gives two components \(u_s, v_s\) of the surface geostrophic current writen in equation (5) The pressure of the level is

\[ \begin{align*} u_s = \frac{g}{f} \frac{\delta \varsigma}{\delta y}\\ v_s = \frac{g}{f} \frac{\delta \varsigma}{\delta x} \tag{5} \end{align*} \]

The Oceanic Topography

We can describe the topography of the sea surface \(\varsigma\) as the height of the sea surface relative to a particular level surface—the geoid. And geoid is the level surface that coincide with the surface of the ocean at rest. Thus, according to equation (eq:eq5) the surface geostrophic currents are proportional to the slope of the topography (figure 2)

cowplot::ggdraw() + cowplot::draw_image("geoid-01.png")

Figure 2: The slope of the sea surface relative to the geoid is direcltyl related to surface geostrophic currents. The slope of 1 meter per 100 kilometers is typical of strong currents

subsetting the Geostrophic Current

In the previous section I narrated a bit of the physics behind the geostrophic current. We now switch to AVISO that collects and process geostrophic currents from satellite. I first chopped the geostrophic current for the tropical Indian Ocean region with griddap() function. Like an R array, ERDDAP grids are subsetted by setting limits on the dimension variables, the difference being that a subset is defined in coordinate space (latitude values, longitude values, time values) rather than array space as is done with R arrays. Thus for erdTAgeo1day_LonPM180 the desired area of the data extract is latitude limits of (35 ^oS, 10 ^oN), longitude limits of (25 ^oE, 65 ^oE), and time limits of (1993-01-01, 2012-12-31) the following would be passed to the function griddap():

#get metadata information on an ERDDAP dataset
geoinfo = info("erdTAgeo1day_LonPM180")

## set spatial extent
lon = c(25,65)
lat =  c(-35,10)
## set temperal extent
time = c("1993-01-01", "2012-12-07")

## download the data
geostrophic = griddap(geoinfo, longitude = lon, latitude = lat, 
                       time = time, fields = c("u_current", "v_current"))

Using equation (6) geostrophic current velocity was computed from the eastward (U) and northward (V) components of sea Water Velocity.

\[ \begin{equation} Velocity (ms^{-1})\:=\: \sqrt{(U^2+V^2)} \tag{6} \end{equation} \]

wio.geostrophic =  geostrophic %>% 
  mutate(date = as.Date(time),
         vel = sqrt(u_current^2 + v_current^2)) %>%
  select(date, lon,lat, u = u_current, v = v_current, vel) %>%
  filter(date >= dmy(010112))

Table 2 show the sample of the geostrophic structure of the zonal (V) , meridional (U), current along the longitude and latitude and time within the tropical Indian Ocean region.

wio %>% sample_n(15) %>%
  kableExtra::kable("html", col.names = c("Date","Longitude", "Latitude", "Zonal (U)", 
                                          "Meridional (V)","Velocity"), 
                    align = "c", digits = 4,   
                    caption = "wind data in the Western Indian Ocean Region")  %>%
  kableExtra::column_spec(column = 1:6, width = "8cm") %>%
  kableExtra::add_header_above(c("", "Geographical Location" = 2, "Current Velocity" = 3))

Climatology geostrophic Current

To calculate the monthly climatology, the date variable in the dataset was separated into three variables—year, month and day using function from lubridate package (Grolemund & Wickham, 2011). Once the month variable were decomposed, the geostrophic values for each grid was grouped into months from January to December. Then the mean geostrophic value in the grid was calculated.

# decompose day, month, and year as separate variable from the date variable
wio = wio %>% 
  mutate(day = day(date), month = month(date) %>%
           as.integer(), year = year(date)%>%
           as.integer()) 
## group current based on longitude, latitude and month
wio.month = wio %>% group_by(lon,lat,month) %>% 
  summarise(u = mean(u, na.rm = T),v = mean(v, na.rm = T), vel= mean(vel, na.rm = T))

## ungroup grouped data frame
wio.month = wio.month %>% ungroup()

Table 2 shows climatological monthly ocean current dynamics in the tropical Indian Ocean region collected between 2012-01-04 and 2012-12-07. Figure 3 show the animated climatological monthly geostrophic current. The current velocity were mapped and animated with gganimate—grammar of animated graphics(Pedersen & Robinson, 2017) and ggplot2—grammar of graphics packages (Wickham, 2016)

ggplot() +
  geom_raster(data = wio.month, 
              aes(x = lon, y = lat, fill = vel), interpolate = TRUE)+ 
  geom_sf(data = spData::world, col = 1, fill = "ivory")+
  coord_sf(xlim = c(30,60), ylim = c(-30,5))+
  # ggrepel::geom_text_repel(data = country.wio, 
  #           aes(x = lon, y = lat, label = name), nudge_x = 0.05)+
  theme_bw() + 
  theme(panel.background = element_rect(fill = "lightblue"),
        axis.text = element_text(size = 11, colour = 1),
        panel.grid = element_line(colour = NA), 
        legend.key.height = unit(3, "lines"),
        legend.position = c(.92,.25),
        legend.background = element_rect(colour = 1))+
  scale_x_continuous(breaks = c(35,55))+
  scale_y_continuous(breaks = c(-25, 2)) +
  labs(x = NULL, y = NULL) +
  scale_fill_gradientn(colours = oceColors9A(120))+
  labs(title = 'Month: {frame_time}') +
  transition_time(month) +
  ease_aes('linear')

Figure 3: Climatology monthly mean Current in the tropical Indian Ocean region

Figure 4 show the animation of monthly climatological geostrophic current in the area around Somali current. We notice the area experience high speed geostrophic current above 0.5 ms^-1 from June to October.

somali = wio.month %>% filter(lon < 53 & lat >2)

ggplot(data =somali , 
       aes(x = month, y = vel, fill = as.factor(month)))+
  geom_boxplot(show.legend = FALSE)+
  scale_fill_brewer(palette = "Paired")+
  cowplot::theme_cowplot()+
  scale_x_continuous(breaks = 1:12, 
                     labels = month(seq(dmy(010118), dmy(311218), by = "month"), 
                                    label = TRUE, abbr = TRUE))  +
  labs(title = "{closest_state}")+
  transition_states(month, transition_length = 5, state_length = 1, wrap = TRUE, states = ) +
  ease_aes('linear') +
  shadow_mark()

Figure 4: The monthly variation of geostrophic current in area around somali current

Remarks

The package rerddap has made access of oceanographic and environmental data into R. Once the data is in the workspace, we can transform, visualize and even animate to uncover subtle clue happening in the tropical Indian Ocean.