diff --git a/ci/requirements/doc.yml b/ci/requirements/doc.yml
index a73c6679322..f722b9d1535 100644
--- a/ci/requirements/doc.yml
+++ b/ci/requirements/doc.yml
@@ -23,6 +23,7 @@ dependencies:
   - pooch
   - pip
   - pydata-sphinx-theme>=0.4.3
+  - pyproj
   - rasterio>=1.1
   - seaborn
   - setuptools
diff --git a/doc/examples/visualization_gallery.ipynb b/doc/examples/visualization_gallery.ipynb
index 49717d53b65..831f162d998 100644
--- a/doc/examples/visualization_gallery.ipynb
+++ b/doc/examples/visualization_gallery.ipynb
@@ -191,7 +191,10 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true,
+    "tags": []
+   },
    "source": [
     "## `imshow()` and rasterio map projections\n",
     "\n",
@@ -213,7 +216,7 @@
     "\n",
     "# The data is in UTM projection. We have to set it manually until\n",
     "# https://github.com/SciTools/cartopy/issues/813 is implemented\n",
-    "crs = ccrs.UTM('18N')\n",
+    "crs = ccrs.UTM('18')\n",
     "\n",
     "# Plot on a map\n",
     "ax = plt.subplot(projection=crs)\n",
@@ -242,20 +245,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from rasterio.warp import transform\n",
+    "from pyproj import Transformer\n",
     "import numpy as np\n",
     "\n",
     "da = xr.tutorial.open_rasterio(\"RGB.byte\")\n",
     "\n",
-    "# Compute the lon/lat coordinates with rasterio.warp.transform\n",
-    "ny, nx = len(da['y']), len(da['x'])\n",
     "x, y = np.meshgrid(da['x'], da['y'])\n",
-    "\n",
-    "# Rasterio works with 1D arrays\n",
-    "lon, lat = transform(da.crs, {'init': 'EPSG:4326'},\n",
-    "                     x.flatten(), y.flatten())\n",
-    "lon = np.asarray(lon).reshape((ny, nx))\n",
-    "lat = np.asarray(lat).reshape((ny, nx))\n",
+    "transformer = Transformer.from_crs(da.crs, \"EPSG:4326\", always_xy=True)\n",
+    "lon, lat = transformer.transform(x, y)\n",
     "da.coords['lon'] = (('y', 'x'), lon)\n",
     "da.coords['lat'] = (('y', 'x'), lat)\n",
     "\n",
@@ -265,14 +262,14 @@
     "# Plot on a map\n",
     "ax = plt.subplot(projection=ccrs.PlateCarree())\n",
     "greyscale.plot(ax=ax, x='lon', y='lat', transform=ccrs.PlateCarree(),\n",
-    "               cmap='Greys_r', add_colorbar=False)\n",
+    "               cmap='Greys_r', shading=\"auto\",add_colorbar=False)\n",
     "ax.coastlines('10m', color='r')"
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -286,7 +283,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.9.7"
   }
  },
  "nbformat": 4,
diff --git a/doc/gallery/plot_rasterio.py b/doc/gallery/plot_rasterio.py
index 99eb1fd1daf..8294e01975f 100644
--- a/doc/gallery/plot_rasterio.py
+++ b/doc/gallery/plot_rasterio.py
@@ -18,7 +18,7 @@
 import cartopy.crs as ccrs
 import matplotlib.pyplot as plt
 import numpy as np
-from rasterio.warp import transform
+from pyproj import Transformer
 
 import xarray as xr
 
@@ -26,14 +26,9 @@
 url = "https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif"
 da = xr.open_rasterio(url)
 
-# Compute the lon/lat coordinates with rasterio.warp.transform
-ny, nx = len(da["y"]), len(da["x"])
-x, y = np.meshgrid(da["x"], da["y"])
-
-# Rasterio works with 1D arrays
-lon, lat = transform(da.crs, {"init": "EPSG:4326"}, x.flatten(), y.flatten())
-lon = np.asarray(lon).reshape((ny, nx))
-lat = np.asarray(lat).reshape((ny, nx))
+# Compute the lon/lat coordinates with pyproj
+transformer = Transformer.from_crs(da.crs, "EPSG:4326", always_xy=True)
+lon, lat = transformer.transform(*np.meshgrid(da["x"], da["y"]))
 da.coords["lon"] = (("y", "x"), lon)
 da.coords["lat"] = (("y", "x"), lat)