diff --git a/.gitignore b/.gitignore
index 31854eb..b9263da 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,3 +12,4 @@ workshop/jupyter/content/data/woudc-stations.xml
 workshop/jupyter/content/data/srtm/srtm_high_pass.geotiff
 workshop/jupyter/content/notebooks/test/10-boreholes.gml
 workshop/jupyter/content/notebooks/test/10-populated-places-ba.json
+workshop/jupyter/content/data/output
diff --git a/workshop/jupyter/content/data/world.rgb.tif.aux.xml b/workshop/jupyter/content/data/world.rgb.tif.aux.xml
new file mode 100644
index 0000000..b4eadeb
--- /dev/null
+++ b/workshop/jupyter/content/data/world.rgb.tif.aux.xml
@@ -0,0 +1,18 @@
+<PAMDataset>
+  <PAMRasterBand band="1">
+    <Description>Red band</Description>
+    <Metadata>
+      <MDI key="STATISTICS_MINIMUM">0</MDI>
+      <MDI key="STATISTICS_MAXIMUM">255</MDI>
+      <MDI key="STATISTICS_MEAN">107.15079498291</MDI>
+      <MDI key="STATISTICS_STDDEV">66.706750840795</MDI>
+      <MDI key="STATISTICS_VALID_PERCENT">100</MDI>
+    </Metadata>
+  </PAMRasterBand>
+  <PAMRasterBand band="2">
+    <Description>Green band</Description>
+  </PAMRasterBand>
+  <PAMRasterBand band="3">
+    <Description>Blue band</Description>
+  </PAMRasterBand>
+</PAMDataset>
diff --git a/workshop/jupyter/content/notebooks/04-vector-data.ipynb b/workshop/jupyter/content/notebooks/04-vector-data.ipynb
index 24bdb77..7a7d040 100644
--- a/workshop/jupyter/content/notebooks/04-vector-data.ipynb
+++ b/workshop/jupyter/content/notebooks/04-vector-data.ipynb
@@ -743,12 +743,14 @@
    },
    "outputs": [],
    "source": [
-    "ds = ogr.Open('../data/04-ogr-out.gml')\n",
-    "layer = ds.GetLayer(0)\n",
-    "print(layer.GetFeatureCount())\n",
-    "print(layer.GetFeature(0).GetField('name'))\n",
-    "print(layer.GetFeature(0).GetField('code'))\n",
-    "ds.Destroy()"
+    "with ogr.Open('../data/04-ogr-out.gml') as ds:\n",
+    "  layer = ds.GetLayer(0)\n",
+    "  ft = layer.GetFeature(0)\n",
+    "  display(ft.GetField('name'))\n",
+    "  display(ft.GetField('code'))\n",
+    "  # print the geometry using shapely\n",
+    "  from shapely import from_wkb\n",
+    "  display(from_wkb(bytes(ft.GetGeometryRef().ExportToWkb())))\n"
    ]
   },
   {
@@ -850,10 +852,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#load it as a GeoDataFrame, i.e. a Pandas DataFrame with with a geometry data column\n",
+    "#load a file as a GeoDataFrame, i.e. a Pandas DataFrame with with a geometry data column\n",
     "countries = gpd.read_file('../data/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp')\n",
     "\n",
-    "countries"
+    "# notice how a geodataframe is printed as a table by jupyter/ipython\n",
+    "countries\n",
+    "# also try `countries.plot()`, geopandas uses matplotlib to print a map of the dataframe"
    ]
   },
   {
diff --git a/workshop/jupyter/content/notebooks/05-raster-data.ipynb b/workshop/jupyter/content/notebooks/05-raster-data.ipynb
index cb2c494..1b5d243 100644
--- a/workshop/jupyter/content/notebooks/05-raster-data.ipynb
+++ b/workshop/jupyter/content/notebooks/05-raster-data.ipynb
@@ -62,11 +62,11 @@
    "source": [
     "from osgeo import gdal\n",
     "gdal.UseExceptions()\n",
-    "ds = gdal.Open('../data/world.rgb.tif')\n",
-    "band = ds.GetRasterBand(1)\n",
-    "print('band:', band.GetDescription())\n",
-    "min_val, max_val, mean_val, std_val = band.ComputeStatistics(False)\n",
-    "print('min:', min_val, 'max:', max_val, 'mean:', mean_val,'std:', std_val)"
+    "with gdal.Open('../data/world.rgb.tif') as ds:\n",
+    "  band = ds.GetRasterBand(1)\n",
+    "  print('band:', band.GetDescription())\n",
+    "  min_val, max_val, mean_val, std_val = band.ComputeStatistics(False)\n",
+    "  print('min:', min_val, 'max:', max_val, 'mean:', mean_val,'std:', std_val)"
    ]
   },
   {
@@ -109,8 +109,8 @@
     "import matplotlib.pyplot as plt\n",
     "# Notice how the zarr file is referenced within the zip file\n",
     "raster_path = '/vsizip/../data/sample_geozarr.zarr.zip/sample_geozarr.zarr'\n",
-    "dataset = gdal.Open(raster_path)\n",
-    "r = plt.imshow(ds.GetRasterBand(1).ReadAsArray())\n"
+    "with gdal.Open(raster_path) as ds:\n",
+    "  r = plt.imshow(ds.GetRasterBand(1).ReadAsArray())\n"
    ]
   },
   {