Merged dev to main for new release

Author: bjorkqvi

dnora/aux_funcs.py

@@ -5,7 +5,9 @@
 from scipy.interpolate import griddata
 from scipy import interpolate
 import os, re, glob
-import pyfimex0 as pyfi
+# This is imported inside the pyfimex function, since we want to avoid to require
+# this dependency if someone wan't to run without wind forcing
+#import pyfimex0 as pyfi
 from typing import TYPE_CHECKING, Tuple, List, Union
 from . import file_module
 if TYPE_CHECKING:
@@ -630,6 +632,7 @@ def create_swan_segment_coords(boundary_mask, lon_edges, lat_edges):
 
 def pyfimex(input_file, output_file, projString, xAxisValues, yAxisValues,
           selectVariables, reduceTime_start, reduceTime_end,ensemble_member=False):
+    import pyfimex0 as pyfi
     r = pyfi.createFileReader('netcdf', input_file)
     inter_ll = pyfi.createInterpolator(r)
     inter_ll.changeProjection(pyfi.InterpolationMethod.BILINEAR,

dnora/bnd/bnd_mod.py

@@ -54,8 +54,8 @@ def import_boundary(self, start_time: str, end_time: str,
         boundary_point_grid = Grid(lon=self.grid.lon_edges(),
                                     lat=self.grid.lat_edges(),
                                     name='boundary_points')
-        boundary_point_grid.set_spacing(nx=self.grid.boundary_nx(),
-                                        ny=self.grid.boundary_ny())
+        #boundary_point_grid.set_spacing(nx=self.grid.boundary_nx(),
+        #                                ny=self.grid.boundary_ny())
         boundary_reader.set_restricted_area(boundary_point_grid)
 
         if write_cache or read_cache:
@@ -64,7 +64,7 @@ def import_boundary(self, start_time: str, end_time: str,
         if read_cache and not cache_empty:
             msg.info('Reading boundary data from cache!!!')
             original_boundary_reader = copy(boundary_reader)
-            boundary_reader = DnoraNc(files=glob.glob(f'{cache_folder}/{cache_name}*'), convention = boundary_reader.convention())
+            boundary_reader = DnoraNc(files=glob.glob(f'{cache_folder}/{cache_name}_bnd.nc'), convention = boundary_reader.convention())
 
         self._history.append(copy(boundary_reader))
 
@@ -88,30 +88,19 @@ def import_boundary(self, start_time: str, end_time: str,
         if read_cache and not cache_empty:
             patch_start, patch_end = aux_funcs.determine_patch_periods(self.time(), start_time, end_time)
             if patch_start:
-                msg.info('Not all data found in cache. Patching from original source...')
-
-                lon_all, lat_all = original_boundary_reader.get_coordinates(start_time)
-                inds = point_picker(self.grid, lon_all, lat_all)
-                bnd_list = [self.data]
-                for t0, t1 in zip(patch_start, patch_end):
-                    time, freq, dirs, spec, lon, lat, source = original_boundary_reader(t0, t1, inds)
-                    bnd_list.append(self.compile_to_xr(time, freq, dirs, spec, lon, lat, source))
-
-                self.data = xr.concat(bnd_list, dim="time").sortby('time')
+                msg.info(
+                    f'Timesteps stored in {cache_name}_bnd.nc does not cover the entire time range. set read_cache to false and rerun.')
+                exit(-1)
 
 
         if write_cache:
             msg.info('Caching data:')
-            for month in self.months():
-                cache_file = f"{cache_name}_{month.strftime('%Y-%m')}.nc"
-                t0 = f"{month.strftime('%Y-%m-01')}"
-                d1 = monthrange(int(month.strftime('%Y')), int(month.strftime('%m')))[1]
-                t1 = f"{month.strftime(f'%Y-%m-{d1}')}"
-                outfile = f'{cache_folder}/{cache_file}'
-                if os.path.exists(outfile):
-                    os.remove(outfile)
-                self.data.sel(time=slice(t0, t1)).to_netcdf(outfile)
-                msg.to_file(outfile)
+            cache_file = f"{cache_name}_bnd.nc"
+            outfile = f'{cache_folder}/{cache_file}'
+            if os.path.exists(outfile):
+                os.remove(outfile)
+            self.data.to_netcdf(outfile)
+            msg.to_file(outfile)
 
         #self.mask = [True]*len(self.x())
 

dnora/bnd/pick.py

@@ -31,8 +31,9 @@ class NearestGridPoint(PointPicker):
     """Choose the nearest grid point to each boundary point in the grid.
     Set a maximum allowed distance using `max_dist` (in km) at instantiation time.
     """
-    def __init__(self, max_dist=None):
+    def __init__(self, max_dist=None, remove_duplicate=False):
         self.max_dist = max_dist
+        self.remove_duplicate = remove_duplicate
         pass
 
     def __call__(self, grid, bnd_lon, bnd_lat):
@@ -51,7 +52,11 @@ def __call__(self, grid, bnd_lon, bnd_lat):
             else:
                 msg.plain('DISCARDED, too far: '+ms)
 
-        inds = np.array(inds)
+        if self.remove_duplicate == True:
+            inds = np.unique(np.array(inds))
+            msg.plain('*** Duplicate spectra are removed ***')
+        else:
+            inds = np.array(inds)
         return inds
 
 class Area(PointPicker):

dnora/bnd/read.py

@@ -108,6 +108,7 @@ def _crop(ds):
             return ds.sel(time=slice(start_time, end_time))
         msg.info(f"Getting boundary spectra from cached netcdf (e.g. {self.files[0]}) from {start_time} to {end_time}")
         ds = xr.open_mfdataset(self.files, preprocess=_crop)
+
         ds = ds.sel(x=inds)
         return ds.time.values, ds.freq.values, ds.dirs.values, ds.spec.values, ds.lon.values, ds.lat.values, ds.source
 

dnora/bnd/read_ec.py

@@ -11,7 +11,7 @@
 import cdsapi
 # Import aux_funcsiliry functions
 from .. import msg
-from ..aux_funcs import create_time_stamps, expand_area, int_list_of_days, int_list_of_months, int_list_of_years
+from ..aux_funcs import create_time_stamps, expand_area, day_list
 
 
 def renormalize_era5_spec(bnd_spec):
@@ -25,53 +25,73 @@ def reshape_bnd_spec(bnd_spec):
     pass
 
 
-def download_era5_from_cds(start_time, end_time, lon, lat, dlon, dlat, folder='dnora_wnd_temp') -> str:
-    """Downloads ERA5 10 m wind data from the Copernicus Climate Data Store for a
+def download_era5_from_cds(start_time, end_time, lon, lat, dlon, dlat, folder='dnora_bnd_temp') -> str:
+    """Downloads ERA5 spectral data from the Copernicus Climate Data Store for a
     given area and time period"""
     start_time = pd.Timestamp(start_time)
     end_time = pd.Timestamp(end_time)
     c = cdsapi.Client()
 
-    filename = f'{folder}/EC_ERA5.nc'
+    days = day_list(start_time, end_time)
+    #tt = pd.date_range(start=start_time,end=end_time)
 
-    years = [f'{y:4.0f}' for y in int_list_of_years(start_time, end_time)]
-    months = [f'{m:02.0f}' for m in int_list_of_months(start_time, end_time)]
-    days = [f'{d:02.0f}' for d in int_list_of_days(start_time, end_time)]
+    filename = f'{folder}/EC_ERA5.nc'
 
     # Create string for dates
-    dates = []
-    for y in years:
-        for m in months:
-            for d in days:
-                dates.append(f'{y}-{m}-{d}')
-    dates = '/'.join(dates)
+    #dates = [days[0].strftime('%Y-%m-%d'), days[-1].strftime('%Y-%m-%d')]
+    #dates = '/'.join(dates)
+    dates = f'{str(start_time)[0:10]}/to/{str(end_time)[0:10]}'
 
 
     cds_command ={
         'class': 'ea',
         'date': dates,
         'direction': '1/2/3/4/5/6/7/8/9/10/11/12/13/14/15/16/17/18/19/20/21/22/23/24',
         'domain': 'g',
+        'area': f'{lat[1]}/{lon[0]}/{lat[0]}/{lon[1]}', # north, west, south, east
+        'grid': f'{dlat}/{dlon}',
         'expver': '1',
         'frequency': '1/2/3/4/5/6/7/8/9/10/11/12/13/14/15/16/17/18/19/20/21/22/23/24/25/26/27/28/29/30',
         'param': '251.140',
         'stream': 'wave',
-        'time': '00:00:00/03:00:00/06:00:00/09:00:00/12:00:00/15:00:00/18:00:00/21:00:00',
-        'area': f'{lat[1]+0.0001}/{lon[0]}/{lat[0]}/{lon[1]+0.0001}', # north, west, south, east
-        'grid': f'{dlon}/{dlat}',
+        'time': '00:00:00/01:00:00/02:00:00/03:00:00/04:00:00/05:00:00/06:00:00/07:00:00/08:00:00/09:00:00/10:00:00/11:00:00/12:00:00/13:00:00/14:00:00/15:00:00/16:00:00/17:00:00/18:00:00/19:00:00/20:00:00/21:00:00/22:00:00/23:00:00',
         'type': 'an',
         'format': 'netcdf',
-        }
+    }
+    #print(cds_command)
+
+    # cds_command ={
+    #     'class': 'ea',
+    #     'date': dates,
+    #     'direction': '/'.join([f'{n+1:01.0f}' for n in range(24)]),
+    #     'domain': 'g',
+    #     'expver': '1',
+    #     'frequency': '/'.join([f'{n+1:01.0f}' for n in range(30)]),
+    #     'param': '251.140',
+    #     'stream': 'wave',
+    #     'time': '00:00:00/03:00:00/06:00:00/09:00:00/12:00:00/15:00:00/18:00:00/21:00:00',
+    #     'area': f'{lat[1]}/{lon[0]}/{lat[0]}/{lon[1]}', # north, west, south, east
+    #     'grid': f'{dlon}/{dlat}',
+    #     'type': 'an',
+    #     'format': 'netcdf',
+    #     }
 
     c.retrieve('reanalysis-era5-complete', cds_command, filename)
     return filename
 class ERA5(BoundaryReader):
+
+    def __init__(self):
+        self.dlon = 0.5
+        self.dlat = 0.5
+
     def convention(self) -> str:
         return 'Ocean'
 
     def get_coordinates(self, start_time) -> Tuple:
         """Reads first time instance of first file to get longitudes and latitudes for the PointPicker"""
-        point_list = self.get_restricted_area()._point_list()
+        restricted_area = self.get_restricted_area()
+        restricted_area.set_spacing(dlon=self.dlon, dlat=self.dlat)
+        point_list = restricted_area._point_list()
         lon_all = point_list[:,0]
         lat_all = point_list[:,1]
 
@@ -87,23 +107,31 @@ def __call__(self, start_time, end_time, inds) -> Tuple:
             os.mkdir(temp_folder)
             print("Creating folder %s..." % temp_folder)
 
-        msg.plain("Removing old files from temporary folder...")
-        for f in glob.glob(f"{temp_folder}/EC_ERA5.nc"):
-            os.remove(f)
+        local_read = False
+
+        if not local_read:
+            msg.plain("Removing old files from temporary folder...")
+            for f in glob.glob(f"{temp_folder}/EC_ERA5.nc"):
+                os.remove(f)
 
         restricted_area = self.get_restricted_area()
+        lon = np.floor(np.array(restricted_area.lon_edges())/self.dlon)*self.dlon
+        lat = np.floor(np.array(restricted_area.lat_edges())/self.dlat)*self.dlat
 
-        nc_file = download_era5_from_cds(start_time, end_time,
-                                        lon=restricted_area.lon_edges(),
-                                        lat=restricted_area.lat_edges(),
-                                        dlon=restricted_area.dlon(),
-                                        dlat=restricted_area.dlat(),
-                                        folder=temp_folder)
 
 
+        if local_read:
+            nc_file = f'{temp_folder}/EC_ERA5.nc'
+        else:
+            nc_file = download_era5_from_cds(start_time, end_time,
+                                            lon=(lon[0]-self.dlon, lon[1]+self.dlon),
+                                            lat=(lat[0]-self.dlat, lat[1]+self.dlat),
+                                            dlon=self.dlon,
+                                            dlat=self.dlat,
+                                            folder=temp_folder)
 
         bnd_spec = xr.open_dataset(nc_file)
-
+        bnd_spec = bnd_spec.sortby("time")
         bnd_spec = renormalize_era5_spec(bnd_spec)
 
         lon, lat = np.meshgrid(bnd_spec.longitude.values, bnd_spec.latitude.values[::-1])