Web Feature Service - Accessing region countours saved on a GeoServer

In this example, we’re going to look at some layers that are currently accessible on our instance of GeoServer. With WFS, we can see what is available, collect the layers we want by using a query, download the results in geoJSON, and visualize them using geopandas.

We begin by loading the libraries needed for parsing and downloading from WFS and for opening and visualizing the results

%matplotlib inline

import os

os.environ["USE_PYGEOS"] = "0"  # force use Shapely with GeoPandas

import geopandas as gpd

# Import WFS from owslib
from owslib.wfs import WebFeatureService
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 os
      5 os.environ["USE_PYGEOS"] = "0"  # force use Shapely with GeoPandas

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:1043, in FontManager.__init__(self, size, weight)
   1040 for path in [*findSystemFonts(paths, fontext=fontext),
   1041              *findSystemFonts(fontext=fontext)]:
   1042     try:
-> 1043         self.addfont(path)
   1044     except OSError as exc:
   1045         _log.info("Failed to open font file %s: %s", path, exc)

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/font_manager.py:1072, in FontManager.addfont(self, path)
   1070 if Path(path).suffix.lower() == ".afm":
   1071     with open(path, "rb") as fh:
-> 1072         font = _afm.AFM(fh)
   1073     prop = afmFontProperty(path, font)
   1074     self.afmlist.append(prop)

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/_afm.py:364, in AFM.__init__(self, fh)
    362 self._header = _parse_header(fh)
    363 self._metrics, self._metrics_by_name = _parse_char_metrics(fh)
--> 364 self._kern, self._composite = _parse_optional(fh)

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/_afm.py:353, in _parse_optional(fh)
    350     key = line.split()[0]
    352     if key in optional:
--> 353         d[key] = optional[key](fh)
    355 return d[b'StartKernData'], d[b'StartComposites']

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/_afm.py:267, in _parse_kern_pairs(fh)
    265     if len(vals) != 4 or vals[0] != b'KPX':
    266         raise RuntimeError('Bad kern pairs line: %s' % line)
--> 267     c1, c2, val = _to_str(vals[1]), _to_str(vals[2]), _to_float(vals[3])
    268     d[(c1, c2)] = val
    269 raise RuntimeError('Bad kern pairs parse')

File ~/checkouts/readthedocs.org/user_builds/pavics-sdi/conda/latest/lib/python3.10/site-packages/matplotlib/_afm.py:69, in _to_str(x)
     68 def _to_str(x):
---> 69     return x.decode('utf8')

KeyboardInterrupt: 

We start by making a connection to the PAVICS instance we have locally on our server. Using WFS, we can very quickly see the contents, which are the layers and the workspaces they’re located with (ie: TravisTest, scratchTJS). These layer names act as dictionaries for

# NBVAL_IGNORE_OUTPUT

wfs_url = "https://pavics.ouranos.ca/geoserver/wfs"  # TEST_USE_PROD_DATA

# Connect to GeoServer WFS service.
wfs = WebFeatureService(wfs_url, version="2.0.0")

# Print the list of available layers
sorted(wfs.contents.keys())
['TravisTest:NE_Admin_Level0',
 'TravisTest:mrc_poly',
 'TravisTest:region_admin_poly',
 'public:CANOPEX_5797_basinBoundaries',
 'public:CANVEC_hydro_waterbodies',
 'public:CanVec_Rivers',
 'public:CanVec_WaterBodies',
 'public:HydroLAKES_points',
 'public:HydroLAKES_poly',
 'public:USGS_HydroBASINS_lake_ar_lev12',
 'public:USGS_HydroBASINS_lake_na_lev12',
 'public:canada_admin_boundaries',
 'public:decamillenial_flood_CC',
 'public:gaspesie_mrc',
 'public:global_admin_boundaries',
 'public:ne_10m_populated_places',
 'public:quebec_admin_boundaries',
 'public:quebec_health_regions',
 'public:quebec_mrc_boundaries',
 'public:quebec_muni_boundaries',
 'public:routing_1kmLakes_07',
 'public:routing_1kmLakes_08',
 'public:routing_1kmLakes_09',
 'public:routing_1kmLakes_10',
 'public:routing_1kmLakes_11',
 'public:routing_1kmLakes_12',
 'public:routing_allLakes_07',
 'public:routing_allLakes_08',
 'public:routing_allLakes_09',
 'public:routing_allLakes_10',
 'public:routing_allLakes_11',
 'public:routing_allLakes_12',
 'public:test_regions',
 'public:test_regions_lambert',
 'public:usa_admin_boundaries',
 'public:wshed_bound_n1',
 'public:wshed_bound_n2',
 'public:wshed_bound_n3']

More information about each layer is stored in the contents dictionary.

# NBVAL_IGNORE_OUTPUT

sorted_layer_ids = list(sorted(wfs.contents.keys()))
canada_admin_boundaries_index = sorted_layer_ids.index("public:canada_admin_boundaries")

for layerID in sorted_layer_ids[
    canada_admin_boundaries_index - 1 : canada_admin_boundaries_index + 2
]:
    layer = wfs[layerID]
    print("Layer ID:", layerID)
    print("Title:", layer.title)
    print("Boundaries:", layer.boundingBoxWGS84, "\n")
Layer ID: public:USGS_HydroBASINS_lake_na_lev12
Title: USGS_HydroBASINS_lake_na_lev12
Boundaries: (-180.0, -90.0, 180.0, 90.0) 

Layer ID: public:canada_admin_boundaries
Title: Canada Administrative Boundaries
Boundaries: (-141.01807315799994, 41.681435425000075, -52.61940850399992, 83.13550252400006) 

Layer ID: public:decamillenial_flood_CC
Title: decamillenial_flood_CC
Boundaries: (-180.0, -90.0, 180.0, 90.0) 

We can then perform a GetFeatures call using the layer name as a target. This returns an IOstream that can be written as a geoJSON file, a common file format for vector data served throughout the web. To reduce the download size, we’ll only fetch the features (here polygons), intersecting a small region defined by a bounding box.

layer_id = "public:canada_admin_boundaries"
meta = wfs.contents[layer_id]
print(meta.title)

# Get the actual data
data = wfs.getfeature(
    typename="public:canada_admin_boundaries",
    bbox=(-74.5, 45.2, -73, 46),
    outputFormat="JSON",
)

# Write to file
fn = "output.geojson"
with open(fn, "wb") as fh:
    fh.write(data.read())
Canada Administrative Boundaries

Once the geoJSON file is downloaded, we can either open it with a GIS application or we can read the features using geopandas.

layers = gpd.read_file(fn)
layers.plot()
<Axes: >
../_images/01125302c7c377dad6ff667ecbc8a30112b41f3230af5745698b22c16e264e94.png