Web Coverage Service - Accessing GeoMet data using owslib

In this notebook we’ll connect to Environment Canada’s GeoMet service and fetch data using the WCS standard.

%matplotlib inline

import matplotlib.pyplot as plt
import xarray as xr
from owslib.wcs import WebCoverageService
Matplotlib is building the font cache; this may take a moment.
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[1], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
      3 import matplotlib.pyplot as plt
      4 import xarray as xr

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/IPython/core/interactiveshell.py:2456, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2454     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2455 with self.builtin_trap:
-> 2456     result = fn(*args, **kwargs)
   2458 # The code below prevents the output from being displayed
   2459 # when using magics with decorator @output_can_be_silenced
   2460 # when the last Python token in the expression is a ';'.
   2461 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/IPython/core/magics/pylab.py:99, in PylabMagics.matplotlib(self, line)
     97     print("Available matplotlib backends: %s" % backends_list)
     98 else:
---> 99     gui, backend = self.shell.enable_matplotlib(args.gui.lower() if isinstance(args.gui, str) else args.gui)
    100     self._show_matplotlib_backend(args.gui, backend)

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3633, in InteractiveShell.enable_matplotlib(self, gui)
   3612 def enable_matplotlib(self, gui=None):
   3613     """Enable interactive matplotlib and inline figure support.
   3614 
   3615     This takes the following steps:
   (...)
   3631         display figures inline.
   3632     """
-> 3633     from matplotlib_inline.backend_inline import configure_inline_support
   3635     from IPython.core import pylabtools as pt
   3636     gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib_inline/__init__.py:1
----> 1 from . import backend_inline, config  # noqa
      2 __version__ = "0.1.6"  # noqa

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib_inline/backend_inline.py:8
      6 import matplotlib
      7 from matplotlib import colors
----> 8 from matplotlib.backends import backend_agg
      9 from matplotlib.backends.backend_agg import FigureCanvasAgg
     10 from matplotlib._pylab_helpers import Gcf

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/backends/backend_agg.py:31
     29 import matplotlib as mpl
     30 from matplotlib import _api, cbook
---> 31 from matplotlib.backend_bases import (
     32     _Backend, FigureCanvasBase, FigureManagerBase, RendererBase)
     33 from matplotlib.font_manager import fontManager as _fontManager, get_font
     34 from matplotlib.ft2font import (LOAD_FORCE_AUTOHINT, LOAD_NO_HINTING,
     35                                 LOAD_DEFAULT, LOAD_NO_AUTOHINT)

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/backend_bases.py:46
     43 import numpy as np
     45 import matplotlib as mpl
---> 46 from matplotlib import (
     47     _api, backend_tools as tools, cbook, colors, _docstring, text,
     48     _tight_bbox, transforms, widgets, is_interactive, rcParams)
     49 from matplotlib._pylab_helpers import Gcf
     50 from matplotlib.backend_managers import ToolManager

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/text.py:16
     14 from . import _api, artist, cbook, _docstring
     15 from .artist import Artist
---> 16 from .font_manager import FontProperties
     17 from .patches import FancyArrowPatch, FancyBboxPatch, Rectangle
     18 from .textpath import TextPath, TextToPath  # noqa # Logically located here

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/font_manager.py:1582
   1578     _log.info("generated new fontManager")
   1579     return fm
-> 1582 fontManager = _load_fontmanager()
   1583 findfont = fontManager.findfont
   1584 get_font_names = fontManager.get_font_names

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/font_manager.py:1576, in _load_fontmanager(try_read_cache)
   1574             _log.debug("Using fontManager instance from %s", fm_path)
   1575             return fm
-> 1576 fm = FontManager()
   1577 json_dump(fm, fm_path)
   1578 _log.info("generated new fontManager")

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/font_manager.py:1041, in FontManager.__init__(self, size, weight)
   1038 try:
   1039     for fontext in ["afm", "ttf"]:
   1040         for path in [*findSystemFonts(paths, fontext=fontext),
-> 1041                      *findSystemFonts(fontext=fontext)]:
   1042             try:
   1043                 self.addfont(path)

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/font_manager.py:280, in findSystemFonts(fontpaths, fontext)
    278     fontpaths = []
    279 else:
--> 280     installed_fonts = _get_fontconfig_fonts()
    281     if sys.platform == 'darwin':
    282         fontpaths = [*X11FontDirectories, *OSXFontDirectories]

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/font_manager.py:258, in _get_fontconfig_fonts()
    255         _log.warning(  # fontconfig 2.7 implemented --format.
    256             'Matplotlib needs fontconfig>=2.7 to query system fonts.')
    257         return []
--> 258     out = subprocess.check_output(['fc-list', '--format=%{file}\\n'])
    259 except (OSError, subprocess.CalledProcessError):
    260     return []

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/subprocess.py:421, in check_output(timeout, *popenargs, **kwargs)
    418         empty = b''
    419     kwargs['input'] = empty
--> 421 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
    422            **kwargs).stdout

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/subprocess.py:505, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
    503 with Popen(*popenargs, **kwargs) as process:
    504     try:
--> 505         stdout, stderr = process.communicate(input, timeout=timeout)
    506     except TimeoutExpired as exc:
    507         process.kill()

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/subprocess.py:1141, in Popen.communicate(self, input, timeout)
   1139     self._stdin_write(input)
   1140 elif self.stdout:
-> 1141     stdout = self.stdout.read()
   1142     self.stdout.close()
   1143 elif self.stderr:

KeyboardInterrupt: 
# NBVAL_IGNORE_OUTPUT
wcs_url = "http://geo.weather.gc.ca/geomet/?lang=en&service=WCS"

# Connect to service
wcs = WebCoverageService(wcs_url, version="2.0.1")
print(wcs.identification.title)

# List some of the content available
sorted(list(wcs.contents.keys()))[:10]
MSC GeoMet — GeoMet-Weather 2.18.5
['GDPS.CONV_KINDEX.PT3H',
 'GDPS.CONV_KINDEX.PT6H',
 'GDPS.CONV_ML-LCL-HGT.3h',
 'GDPS.CONV_ML-LCL-HGT.6h',
 'GDPS.CONV_ML-LI.400.3h',
 'GDPS.CONV_ML-LI.400.6h',
 'GDPS.CONV_ML-LI.500.3h',
 'GDPS.CONV_ML-LI.500.6h',
 'GDPS.CONV_ML-LI.600.3h',
 'GDPS.CONV_ML-LI.600.6h']

Now let’s get some information about a given layer, here the salinity.

layerid = "OCEAN.GIOPS.3D_SALW_0000"
temp = wcs[layerid]
# Title
print("Layer title :", temp.title)
# bounding box
print("BoundingBox :", temp.boundingBoxWGS84)
# supported data formats - we'll use geotiff
print("Formats :", temp.supportedFormats)
# grid dimensions
print("Grid upper limits :", temp.grid.highlimits)
Layer title : None
BoundingBox : None
Formats : ['image/tiff', 'image/x-aaigrid', 'image/netcdf', 'image/png', 'image/jpeg', 'image/png; mode=8bit', 'image/vnd.jpeg-png', 'image/vnd.jpeg-png8']
Grid upper limits : ['1799', '849']

To request data, we need to call the getCoverage service, which requires us specifying the geographic projection, the bounding box, the resolution and format of the output. With GeoMet 2.0.1, we can now get layers in the netCDF format.

format_wcs = "image/netcdf"
bbox_wcs = temp.boundingboxes[0]["bbox"]  # Get the entire domain
crs_wcs = temp.boundingboxes[0]["nativeSrs"]  # Coordinate system
w = int(temp.grid.highlimits[0])
h = int(temp.grid.highlimits[1])

print("Format:", format_wcs)
print("Bounding box:", bbox_wcs)
print("Projection:", crs_wcs)
print(f"Resolution: {w} x {h}")

output = wcs.getCoverage(
    identifier=[
        layerid,
    ],
    crs=crs_wcs,
    bbox=bbox_wcs,
    width=w,
    height=h,
    format=format_wcs,
)
Format: image/netcdf
Bounding box: (-80.1, -180.0, 89.9, 180.0)
Projection: http://www.opengis.net/def/crs/EPSG/0/4326
Resolution: 1799 x 849

We then save the output to disk, open it normally using xarray and plot it’s variable.

fn = layerid + ".nc"
with open(fn, "wb") as fh:
    fh.write(output.read())
ds = xr.open_dataset(fn)
print(ds.data_vars)
ds.Band1.plot()
plt.show()
Data variables:
    Band1    (lat, lon) float32 ...
../_images/a4828a5394aee28735ee026886bcfb651a43ad5078be0d44ff45a664c43b7fcb.png