I am looking into the CF-standards for profile data to implement a netCDF save option for python-ctd.

Until now I was using a customized HDF5 file that, even though it was practical, the file was useless for software like iris. Iris obeys the CF-standards rules and rejects data that does not comply to it. Here is an iris example for plotting profile data.

I decided to try it with a netCDF that was created with the CF profile standards to see how easy it would be to load and plot it once I implement the rules in python-ctd.

Here is the result:

First let's load sea_water_sality cube from a profile example file provided by NODC.

In [2]:

import iris

url = 'http://data.nodc.noaa.gov/thredds/dodsC/testdata/netCDFTemplateExamples/profile/wodStandardLevels.nc'
sal = iris.load_cube(url, 'sea_water_salinity')

print(sal)

sea_water_salinity / (unknown)      (-- : 169; altitude: 26)
     Dimension coordinates:
          altitude                      -              x
     Auxiliary coordinates:
          latitude                      x              -
          longitude                     x              -
          time                          x              -
     Attributes:
          Conventions: CF-1.6
          DODS.dimName: strnlen
          DODS.strlen: 170
          cdm_data_type: Profile
          featureType: profile
          standard_name_vocabulary: CF-1.6

The dimension coordinate has 169 profiles by 26 altitude (positive down, AKA depth!), and the auxiliary coordinates are in the 169 Dimension coordinate. That means I have 169 profiles and the positions are stored in the auxiliary coordinates. Very different from model data were the Dimension coordinate are usually the positions.

Before we try to extract a profile let's plot the data positions to understand what the profile rules are all about.

In [3]:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

from cartopy.feature import LAND, COASTLINE
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

x = sal.coord('longitude').points
y = sal.coord('latitude').points

fig, ax = plt.subplots(figsize=(11, 13),
                       subplot_kw=dict(projection=ccrs.PlateCarree()))
ax.set_extent([-180, 180, -90, 90])
ax.stock_img()
ax.plot(x, y, 'go')

ax.add_feature(COASTLINE)
kw = dict(linewidth=1.5, color='gray', alpha=0.5, linestyle='--')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, **kw)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER

from shapely.geometry.polygon import LinearRing
lons = [105, 105, 125, 125]
lats = [-40, -20, -20, -40]
ring = LinearRing(list(zip(lons, lats)))
sq = ax.add_geometries([ring], ccrs.PlateCarree(), facecolor='none', edgecolor='red')

OK, so we want that one profile inside the square! Let's build the iris Constraint for that box. We will also limit the depth of the of the data to 30 m avoid loading missing data from below it. (That is a shallow profile, but it serve our purpose.)

In [4]:

lon = iris.Constraint(longitude=lambda l:min(lons) < l < max(lons))
lat = iris.Constraint(latitude=lambda l:min(lats) < l < max(lats))
alt = iris.Constraint(altitude=lambda a: a <= 30)

sal_profile = sal.extract(alt & lon & lat)

print(sal_profile)

sea_water_salinity / (unknown)      (altitude: 4)
     Dimension coordinates:
          altitude                           x
     Scalar coordinates:
          latitude: -28.7667 degrees
          longitude: 114.383 degrees
          time: 1980-03-10 00:00:00
     Attributes:
          Conventions: CF-1.6
          DODS.dimName: strnlen
          DODS.strlen: 170
          cdm_data_type: Profile
          featureType: profile
          standard_name_vocabulary: CF-1.6

And finally the profile plot:

In [5]:

import iris.plot as iplt

lon = sal_profile.coord(axis='X').points.squeeze()
lat = sal_profile.coord(axis='Y').points.squeeze()
depth = sal_profile.coord(axis='Z').points.max()

fig, ax = plt.subplots(figsize=(5, 6))
kw = dict(linewidth=2,  color='k', marker='o')
iplt.plot(sal_profile, sal_profile.coord('altitude'), **kw)
ax.grid()
ax.set_ylabel('Depth (m)')
ax.set_xlabel('Salinity (unknown)')
t = ax.set_title('lon: %s\nlat: %s\nMax depth = %s' % (lon, lat, depth))

Soon python-ctd will be saving files according the CF profile standards, and we will be using iris to expore the data.

In [6]:

HTML(html)

Out[6]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

CF-standards for oceanographic profile data

Comments