diff --git a/pyproject.toml b/pyproject.toml
index 1b32946..51ea922 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,6 +34,7 @@ dependencies = [
     "matplotlib >= 3.8.0",
     "netCDF4 >= 1.7",
     "freeqdsk >= 0.4.0",
+    "shapely",
 ]
 dynamic = ["version"]
 
diff --git a/requirements.txt b/requirements.txt
index 0888706..991a969 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,2 +1,2 @@
-# the following line installs requirements from setup.py
+# the following line installs requirements from pyproject.toml
 -e .
diff --git a/zoidberg/poloidal_grid.py b/zoidberg/poloidal_grid.py
index 290e91b..3483f7f 100644
--- a/zoidberg/poloidal_grid.py
+++ b/zoidberg/poloidal_grid.py
@@ -21,6 +21,7 @@
 import warnings
 
 import numpy as np
+import shapely
 from numpy import linspace, pi, zeros
 from scipy.interpolate import RectBivariateSpline
 from scipy.spatial import cKDTree as KDTree
@@ -281,17 +282,26 @@ def __init__(self, R, Z):
         self.R = R
         self.Z = Z
 
-        # Create a KDTree for quick lookup of nearest points
-        n = R.size
-        data = np.concatenate((R.reshape((n, 1)), Z.reshape((n, 1))), axis=1)
-        self.tree = KDTree(data)
-
-        # Create splines for quick interpolation of coordinates
         nx, nz = R.shape
-
         self.nx = nx
         self.nz = nz
 
+        if nx > 4:
+            inner = shapely.Polygon(zip(R[2, :], Z[2, :]))
+            outer = shapely.Polygon(zip(R[-2, :], Z[-2, :]))
+
+            assert (
+                inner.area <= outer.area
+            ), f"""You are trying to create a grid with inner boundary at high x
+and outer boundary at low x. This is against the convention -
+switch inner and outer boundary.
+                {inner.area} {outer.area}"""
+
+        # Create a KDTree for quick lookup of nearest points
+        data = np.concatenate((R.reshape((-1, 1)), Z.reshape((-1, 1))), axis=1)
+        self.tree = KDTree(data)
+
+        # Create splines for quick interpolation of coordinates
         xinds = np.arange(nx)
         zinds = np.arange(nz * 3)
         # Repeat the data in z, to approximate periodicity