automatically convert non lonlat to lonlat

ecoval/gridded.py

@@ -8,7 +8,7 @@
 import nctoolkit as nc
 
 from ecoval.fixers import tidy_warnings
-from ecoval.utils import extension_of_directory, get_extent
+from ecoval.utils import extension_of_directory, get_extent, is_latlon, get_resolution
 from ecoval.session import session_info
 
 
@@ -709,6 +709,22 @@ def gridded_matchup(
                     ds_surface.top()
                 if lon_lim is not None and lat_lim is not None:
                     ds_surface.subset(lon=lon_lim, lat=lat_lim)
+
+
+                ds_surface.run()
+
+                regrid_later = False
+                if is_latlon(ds_surface[0]) is False:
+                    extent = get_extent(ds_surface[0])
+                    lons = [extent[0], extent[1]]
+                    lats = [extent[2], extent[3]]
+                    resolution = get_resolution(ds_surface[0])
+                    lon_res = resolution[0]
+                    lat_res = resolution[1]
+                    ds_surface.to_latlon(lon=lons, lat=lats, res=[lon_res, lat_res], method = "nn")
+                    regrid_later = True
+
+
                 ds_surface.to_nc(out_file, zip=True, overwrite=True)
 
                 # now do the masking etc.
@@ -734,6 +750,9 @@ def gridded_matchup(
                     if lon_lim is not None and lat_lim is not None:
                         ds_obs_annual.subset(lon=lon_lim, lat=lat_lim)
 
+                    if regrid_later:
+                        ds_obs_annual.to_latlon(lon=lons, lat=lats, res=[lon_res, lat_res], method = "nn")
+
                     ds_obs_annual.to_nc(out_file, zip=True, overwrite=True)
 
 

ecoval/utils.py

@@ -1,6 +1,8 @@
 import os
 import nctoolkit as nc
+import xarray as xr
 import numpy as np
+import subprocess
 
 session = dict()
 
@@ -124,3 +126,41 @@ def extension_of_directory(starting_directory, exclude=[]):
     for i in range(levels):
         new_directory = new_directory + "/**"
     return new_directory + "/"
+
+
+
+def is_latlon(ff):
+    ds = nc.open_data(ff, checks = False)
+
+    cdo_result = subprocess.run(
+            f"cdo griddes {ff}",
+            shell=True,
+            stdout=subprocess.PIPE,
+            stderr=subprocess.PIPE,
+        )
+    return "lonlat" in cdo_result.stdout.decode("utf-8")
+
+def get_resolution(ff):
+    ds = nc.open_data(ff, checks = False)
+    var = ds.variables[0]
+    ds = xr.open_dataset(ff)
+    lon_name = [x for x in list(ds.coords) if 'lon' in x][0]
+    lat_name = [x for x in list(ds.coords) if 'lat' in x][0]
+    var_dims =  ds[var].dims
+    extent = get_extent(ff)
+    if lon_name in var_dims:
+        if lat_name in var_dims:
+            n_lon = len(ds[lon_name])
+            n_lat = len(ds[lat_name])
+            lon_res = (extent[1] - extent[0]) / n_lon
+            lat_res = (extent[3] - extent[2]) / n_lat
+            return [lon_res, lat_res]
+    else:
+        # get the final two var_dims
+        var_dims = var_dims[-2:]
+        # lat should be the second
+        n_lat = len(ds[var_dims[1]])
+        n_lon = len(ds[var_dims[0]])
+        lon_res = (extent[1] - extent[0]) / n_lon
+        lat_res = (extent[3] - extent[2]) / n_lat
+        return [lon_res, lat_res]