From 514b560adf7b9094081408da49ecb3e19ef0353b Mon Sep 17 00:00:00 2001
From: Nic Hannah <nichannah@pm.me>
Date: Fri, 1 Apr 2022 11:46:51 +1100
Subject: [PATCH] ERA5 grid support.

---
 esmgrids/base_grid.py |  9 +++++++++
 esmgrids/era5_grid.py | 28 ++++++++++++++++++++++++++++
 esmgrids/util.py      |  2 ++
 3 files changed, 39 insertions(+)
 create mode 100644 esmgrids/era5_grid.py

diff --git a/esmgrids/base_grid.py b/esmgrids/base_grid.py
index cdc6f1f..2a30930 100644
--- a/esmgrids/base_grid.py
+++ b/esmgrids/base_grid.py
@@ -98,12 +98,21 @@ def __init__(self, **kwargs):
         # pole to be fixed.
         self.fix_pole_holes()
 
+        # The base class may implement this if the corner latitudes overrun the
+        # poles
+        self.fix_pole_overruns()
+
         if self.area_t is None and self.calc_areas:
             self.area_t = calc_area_of_polygons(self.clon_t, self.clat_t)
 
     def fix_pole_holes(self):
         pass
 
+
+    def fix_pole_overruns(self):
+        pass
+
+
     @classmethod
     def fromgrid(cls, grid):
         """
diff --git a/esmgrids/era5_grid.py b/esmgrids/era5_grid.py
new file mode 100644
index 0000000..d411a57
--- /dev/null
+++ b/esmgrids/era5_grid.py
@@ -0,0 +1,28 @@
+
+import numpy as np
+import netCDF4 as nc
+
+from .base_grid import BaseGrid
+
+
+class Era5Grid(BaseGrid):
+
+    def __init__(self, h_grid_def, description='ERA5 regular grid'):
+        self.type = 'Arakawa A'
+        self.full_name = 'ERA5'
+
+        with nc.Dataset(h_grid_def) as f:
+            x_t = f.variables['longitude'][:]
+            # ERA5 puts latitudes the wrong way around
+            y_t = np.flipud(f.variables['latitude'][:])
+
+        super(Era5Grid, self).__init__(x_t=x_t, y_t=y_t, calc_areas=False,
+                                        description=description)
+
+    def fix_pole_overruns(self):
+        """
+        ERA5 has grid cell centres at latitudes of -90, 90.
+        """
+
+        self.clat_t[np.where(self.clat_t > 90)] = 90
+        self.clat_t[np.where(self.clat_t < -90)] = -90
diff --git a/esmgrids/util.py b/esmgrids/util.py
index f7a7ad3..5b4e82d 100644
--- a/esmgrids/util.py
+++ b/esmgrids/util.py
@@ -14,6 +14,8 @@ def calc_area_of_polygons(clons, clats):
     and then calculate the area of flat polygon.
 
     This is slow we should do some caching to avoid recomputing.
+
+    FIXME: compare against a Haversine method for speed and accuracy
     """
 
     areas = np.zeros(clons.shape[1:])