diff --git a/interpies/colors.py b/interpies/colors.py
index ff2e3bb..e98a58b 100644
--- a/interpies/colors.py
+++ b/interpies/colors.py
@@ -10,7 +10,7 @@
         Originally recovered from: https://eisbox.net/blog/2011/12/16/windows-8-bit-256-color-palette/
 
 @author: Joseph Barraud
-Geophysics Labs, 2017
+Geophysics Labs, 2017-2024
 """
 
 from matplotlib.colors import LinearSegmentedColormap
diff --git a/interpies/spatial.py b/interpies/spatial.py
index 33fd198..4ff62b7 100644
--- a/interpies/spatial.py
+++ b/interpies/spatial.py
@@ -5,19 +5,21 @@
     Functions to manipulate spatial data (grids and points)
 
 @author: Joseph Barraud
-Geophysics Labs, 2017
+Geophysics Labs, 2017-2024
 """
 
 import subprocess
 
-# import numpy
 import numpy as np
+from numpy.typing import ArrayLike
 
 # import GDAL modules
-from osgeo import ogr, osr
+from osgeo import gdal, ogr, osr
 
+gdal.UseExceptions()  # Enable exceptions
 
-def project_points(inputPoints, s_srs=4326, t_srs=23029):
+
+def project_points(inputPoints: ArrayLike, s_srs: int = 4326, t_srs: int = 23029):
     """
     Reproject a set of points from one spatial reference to another.
 
@@ -46,8 +48,8 @@ def project_points(inputPoints, s_srs=4326, t_srs=23029):
 
     # Loop through the points
     outputPoints = []
-    for XY in inputPoints:
-        point = ogr.CreateGeometryFromWkt('POINT ({} {})'.format(*XY))
+    for xy in inputPoints:
+        point = ogr.CreateGeometryFromWkt('POINT ({} {})'.format(*xy))
         point.Transform(coordTrans)
         outputPoints.append([point.GetX(), point.GetY()])
 
diff --git a/pyproject.toml b/pyproject.toml
index a590f9a..f645b9c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -25,7 +25,7 @@ dependencies = [
     "matplotlib",
     "numpy",
     "gdal==3.8.4",
-    "rasterio>=1.1",
+    "rasterio>=1.3,<1.4",
     "scikit-image",
     "scikit-learn",
     "scipy",
diff --git a/tests/test_spatial.py b/tests/test_spatial.py
new file mode 100644
index 0000000..e5c079a
--- /dev/null
+++ b/tests/test_spatial.py
@@ -0,0 +1,23 @@
+"""Tests for the spatial module."""
+
+from numpy.testing import assert_allclose
+
+from interpies.spatial import project_points
+
+
+def test_project_points():
+    projected = project_points([[45.0, -6], [45.0, -9], [45.0, -12.0], [60, 0.0]], s_srs=4326, t_srs=23029)
+    desired = [
+        [736557.69131476, 4987548.4788836],
+        [500108.90395927, 4983169.23651159],
+        [263542.99025614, 4987422.43017289],
+        [1001073.0017133, 6685825.55236939],
+    ]
+    assert_allclose(projected, desired)
+
+    unprojected = project_points(
+        desired,
+        s_srs=23029,
+        t_srs=4326,
+    )
+    assert_allclose(unprojected, [[45.0, -6], [45.0, -9], [45.0, -12.0], [60, 0.0]], atol=1e-6)