From 04d34531298c818c12a5a421bf058b9c43386d49 Mon Sep 17 00:00:00 2001
From: James Clark <jcl@pml.ac.uk>
Date: Thu, 14 Dec 2023 10:09:57 +0000
Subject: [PATCH] Enable saving of FVCOM element areas

---
 pylag/grid_metrics.pyx | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/pylag/grid_metrics.pyx b/pylag/grid_metrics.pyx
index 6f50e57..c405460 100644
--- a/pylag/grid_metrics.pyx
+++ b/pylag/grid_metrics.pyx
@@ -612,6 +612,23 @@ def create_fvcom_grid_metrics_file(fvcom_file_name: str,
     nbe_data = nbe_data.T
     nv_data = nv_data.T
     sort_adjacency_array(nv_data, nbe_data)
+
+    # Save element areas
+    print('\nCalculating element areas ', end='... ')
+    lon_nodes = fvcom_dataset.variables[lon_var_name][:].astype(DTYPE_FLOAT)
+    lat_nodes = fvcom_dataset.variables[lat_var_name][:].astype(DTYPE_FLOAT)
+    areas = compute_element_areas(nv_data, lon_nodes, lat_nodes,
+                                  coordinate_system='geographic')
+    area_attrs = {'standard_name' : 'areas',
+                  'units' : 'm^2',
+                  'long_name' : 'Element areas'}
+ 
+    # Save areas
+    gm_file_creator.create_variable('area', areas, ('element',),
+            DTYPE_FLOAT, attrs=area_attrs)
+    print('done')
+
+    # Transpose back to original shape as expected by PyLag
     nbe_data = nbe_data.T
     nv_data = nv_data.T