diff --git a/docs/projections.py b/docs/projections.py
index 35fef280f..fb71a5aa8 100644
--- a/docs/projections.py
+++ b/docs/projections.py
@@ -215,23 +215,24 @@
 # Plotting geographic data
 # ------------------------
 #
-# In ProPlot, plotting in `~proplot.axes.GeoAxes` looks pretty much exactly the
-# same as plotting in `~proplot.axes.CartesianAxes`. While cartopy and basemap
-# assume your data is in "map projection" coordinates unless specified otherwise,
-# ProPlot makes longitude-latitude (i.e. Plate Carrée) coordinates the *default*
-# coordinate system for your datasets by passing ``transform=ccrs.PlateCarree()``
-# to cartopy plotting methods and ``latlon=True`` to basemap plotting methods.
-#
-# There are also a couple plotting features specific to `~proplot.axes.GeoAxes`.
-# To ensure a 2D plot like `~matplotlib.axes.Axes.contour` covers the entire globe,
+# In ProPlot, plotting in `~proplot.axes.GeoAxes` is not much different from
+# plotting in `~proplot.axes.CartesianAxes`. ProPlot makes longitude-latitude
+# (i.e. Plate Carrée) coordinates the *default* coordinate system for your
+# datasets by passing ``transform=ccrs.PlateCarree()`` to cartopy plotting
+# commands and ``latlon=True`` to basemap plotting commands. And again, basemap
+# plotting commands are invoked from the `proplot.axes.Axes.GeoAxes` rather
+# than from the `~mpl_toolkits.basemap.Basemap` instance.
+#
+# To ensure 2D plots like `~matplotlib.axes.Axes.contour` cover the entire globe,
 # pass ``globe=True`` to the plotting command. This interpolates your data
-# to the poles and the longitude seams before plotting.
-#
-# To mask out the parts of your data over the land or the ocean, toggle
-# the :rcraw:`land` or and :rcraw:`ocean` settings and make sure the corresponding
-# `zorder <https://matplotlib.org/3.1.1/gallery/misc/zorder_demo.html>`__
-# is high enough to sit above all plotted content, e.g. with
-# ``plot.rc.update({'land': True, 'land.zorder': 5})``.
+# to the poles and across the longitude seams before plotting, having the same
+# effect as cartopy's `~cartopy.util.add_cyclic_point` function and basemap's
+# `~mpl_toolkits.basemap.addcyclic` function.
+#
+# Geographic feature can be drawn underneath data or on top of data by changing the
+# corresponding `zorder <https://matplotlib.org/3.1.1/gallery/misc/zorder_demo.html>`__
+# rc setting. For example, to draw land patches on top of all plotted content as
+# a "land mask," use ``ax.format(land=True, landzorder=4)``.
 # See the :ref:`next section <ug_geoformat>` for details.
 
 # %%
@@ -248,10 +249,18 @@
 
 # Plot data both without and with globe=True
 for globe in (False, True,):
+    string = 'with' if globe else 'without'
     fig, axs = plot.subplots(
         ncols=2, nrows=2, axwidth=2.5,
         proj='kav7', basemap={(1, 3): False, (2, 4): True}
     )
+    axs.format(
+        suptitle=f'Geophysical data {string} global coverage',
+        collabels=['Cartopy example', 'Basemap example'],
+        rowlabels=['Contourf', 'Pcolor'],
+        abc=True, abcstyle='a)', abcloc='ul', abcborder=False,
+        land=True, lonlines=90,
+    )
     for i, ax in enumerate(axs):
         cmap = ('sunset', 'sunrise')[i % 2]
         if i < 2:
@@ -259,14 +268,6 @@
             fig.colorbar(m, loc='b', span=i + 1, label='values', extendsize='1.7em')
         else:
             ax.pcolor(lon, lat, data, cmap=cmap, globe=globe, extend='both')
-    string = 'with' if globe else 'without'
-    axs.format(
-        suptitle=f'Geophysical data {string} global coverage',
-        collabels=['Cartopy example', 'Basemap example'],
-        rowlabels=['Contourf', 'Pcolor'],
-        abc=True, abcstyle='a)', abcloc='ul', abcborder=False,
-        land=True, landzorder=3, lonlines=90,
-    )
 
 
 # %% [raw] raw_mimetype="text/restructuredtext"
diff --git a/proplot/axes/geo.py b/proplot/axes/geo.py
index cfe57dedf..929fcfc7a 100644
--- a/proplot/axes/geo.py
+++ b/proplot/axes/geo.py
@@ -941,10 +941,11 @@ def _update_features(self):
                 + ', '.join(map(repr, constructor.CARTOPY_RESOS)) + '.'
             )
         for name, args in constructor.CARTOPY_FEATURES.items():
+            # Draw feature or toggle feature off
             b = rc.get(name, context=True)
             attr = f'_{name}_feature'
+            feat = getattr(self, attr, None)
             if b is not None:
-                feat = getattr(self, attr, None)
                 if not b:
                     if feat is not None:  # toggle existing feature off
                         feat.set_visible(False)
@@ -953,17 +954,23 @@ def _update_features(self):
                     if not drawn:
                         feat = cfeature.NaturalEarthFeature(*args, reso)
                         feat = self.add_feature(feat)  # convert to FeatureArtist
-                    # For 'lines', need to specify edgecolor and facecolor
-                    # See: https://github.com/SciTools/cartopy/issues/803
-                    kw = rc.category(name, context=drawn)
-                    if name in ('coast', 'rivers', 'borders', 'innerborders'):
-                        kw.update({'edgecolor': kw.pop('color'), 'facecolor': 'none'})
-                    else:
-                        kw.update({'linewidth': 0})
-                    # Update artist attributes (_kwargs used back to v0.5)
-                    # feat.update(kw)  # TODO: check this fails
+
+            # Update artist attributes (FeatureArtist._kwargs used back to v0.5).
+            # For 'lines', need to specify edgecolor and facecolor
+            # See: https://github.com/SciTools/cartopy/issues/803
+            if feat is not None:
+                kw = rc.category(name, context=drawn)
+                if name in ('coast', 'rivers', 'borders', 'innerborders'):
+                    kw.update({'edgecolor': kw.pop('color'), 'facecolor': 'none'})
+                else:
+                    kw.update({'linewidth': 0})
+                if 'zorder' in kw:
+                    # NOTE: Necessary to update zorder directly because _kwargs
+                    # attributes are not applied until draw()... at which point
+                    # matplotlib is drawing in the order based on the *old* zorder.
+                    feat.set_zorder(kw['zorder'])
+                if hasattr(feat, '_kwargs'):
                     feat._kwargs.update(kw)
-                    setattr(self, attr, feat)
 
     def _update_gridlines(
         self, gl, which='major', longrid=None, latgrid=None, nsteps=None,
@@ -1347,8 +1354,8 @@ def _update_features(self):
         for name, method in constructor.BASEMAP_FEATURES.items():
             b = rc.get(name, context=True)
             attr = f'_{name}_feature'
+            feat = getattr(self, attr, None)
             if b is not None:
-                feat = getattr(self, attr, None)
                 if not b:
                     if feat is not None:  # toggle existing feature off
                         for obj in feat:
@@ -1357,12 +1364,13 @@ def _update_features(self):
                     drawn = feat is not None  # if exists, apply *updated* settings
                     if not drawn:
                         feat = getattr(self.projection, method)(ax=self)
-                    kw = rc.category(name, context=drawn)
                     if not isinstance(feat, (list, tuple)):  # list of artists?
                         feat = (feat,)
-                    for obj in feat:
-                        obj.update(kw)
                     setattr(self, attr, feat)
+            if feat is not None:
+                kw = rc.category(name, context=drawn)
+                for obj in feat:
+                    obj.update(kw)
 
     def _update_gridlines(
         self, which='major', longrid=None, latgrid=None, lonarray=None, latarray=None,