Merge pull request #17 from MET-OM/dev_norkyst

Add norkyst800

docs/index.rst

@@ -96,6 +96,9 @@ Several options for **product** are available. Please check the data catalog for
 * For coastal wave NORA3 data: product='NORAC_wave'
 
   Dataset: https://thredds.met.no/thredds/catalog/norac_wave/field/catalog.html
+* For ocean data (sea level, temperature, currents, salinity over depth ) Norkyst800 data (from 2016-09-14 to today): product='NORKYST800'
+
+  Dataset: https://thredds.met.no/thredds/fou-hi/norkyst800v2.html
 
 * For global reananalysis ERA5 (wind and waves): product='ERA5' 
 

examples/example_import_NORA3.py

@@ -8,12 +8,12 @@
 #df_ts = ts.TimeSeries(lon=1.320, lat=53.324,start_time='2000-01-01', end_time='2001-03-31' , product='NORA3_stormsurge')
 #df_ts = ts.TimeSeries(lon=1.320, lat=53.324,start_time='2021-01-01', end_time='2021-03-31' , product='NORA3_atm_sub')
 #df_ts = ts.TimeSeries(lon=3.7, lat=61.8, start_time='2023-01-01', end_time='2023-02-01', product='NORA3_atm3hr_sub')
+#df_ts = ts.TimeSeries(lon=3.73, lat=64.60,start_time='2019-02-26', end_time='2019-02-27' , product='NORKYST800')
 
 
 # Import data from thredds.met.no and save it as csv
 df_ts.import_data(save_csv=True)
 
-#print(df_ts.data)
 # Load data from a local csv-file
 #df_ts.load_data(local_file=df_ts.datafile)
 

metocean_api/ts/aux_funcs.py

@@ -32,8 +32,10 @@ def find_nearest_rotCoord(lon_model, lat_model, lon0, lat0):
 def find_nearest_cartCoord(lon_model, lat_model, lon0, lat0):
     #print('find nearest point...')
     dx = distance_2points(lat0, lon0, lat_model, lon_model)
-    y0 = dx.where(dx == dx.min(), drop=True).y
-    x0 = dx.where(dx == dx.min(), drop=True).x
+    x_name = 'x' if 'x' in dx.dims else 'X' if 'X' in dx.dims else None
+    x0 = dx.where(dx == dx.min(), drop=True)[x_name]    
+    y_name = 'y' if 'y' in dx.dims else 'Y' if 'Y' in dx.dims else None
+    y0 = dx.where(dx == dx.min(), drop=True)[y_name]    
     return x0, y0
 
 def find_nearest_rhoCoord(lon_model, lat_model, lon0, lat0):
@@ -74,6 +76,13 @@ def get_url_info(product, date):
         #infile = 'https://thredds.met.no/thredds/dodsC/nora3_subset_sealevel/sealevel/zeta_nora3era5_N4_'+date.strftime('%Y')+'.nc'
         x_coor_str = 'eta_rho'
         y_coor_str = 'xi_rho'
+    elif product == 'NORKYST800':     
+        if date>=pd.Timestamp('2016-09-14 00:00:00') and date<=pd.Timestamp('2019-02-26 00:00:00'):
+            infile = 'https://thredds.met.no/thredds/dodsC/sea/norkyst800mv0_1h/NorKyst-800m_ZDEPTHS_his.an.'+date.strftime('%Y%m%d%H')+'.nc'
+        elif date>pd.Timestamp('2019-02-26 00:00:00'):
+            infile = 'https://thredds.met.no/thredds/dodsC/fou-hi/norkyst800m-1h/NorKyst-800m_ZDEPTHS_his.an.'+date.strftime('%Y%m%d%H')+'.nc'
+        x_coor_str = 'X'
+        y_coor_str = 'Y'
     elif product.startswith('E39'):     
         infile = 'https://thredds.met.no/thredds/dodsC/obs/buoy-svv-e39/'+date.strftime('%Y/%m/%Y%m')+'_'+product+'.nc'
         x_coor_str = 'longitude'
@@ -92,6 +101,8 @@ def get_date_list(product, start_date, end_date):
         date_list = pd.date_range(start=start_date , end=end_date, freq='MS')
     elif product == 'NORA3_stormsurge':
         date_list = pd.date_range(start=start_date , end=end_date, freq='YS')
+    elif product == 'NORKYST800':
+        date_list = pd.date_range(start=start_date , end=end_date, freq='D')
     elif product.startswith('E39'):
         date_list = pd.date_range(start=start_date , end=end_date, freq='MS')
     return date_list
@@ -113,6 +124,8 @@ def drop_variables(product):
         drop_var = ['longitude','latitude']  
     elif product == 'NORA3_stormsurge':
         drop_var = ['lon_rho','lat_rho']  
+    elif product == 'NORKYST800':
+        drop_var = ['lon','lat']  
     elif product.startswith('E39'):
         drop_var = ['longitude','latitude']
     return drop_var
@@ -144,6 +157,12 @@ def get_near_coord(infile, lon, lat, product):
         lat_near = ds.lat_rho.sel(eta_rho=eta_rho, xi_rho=xi_rho).values[0][0]
         x_coor = eta_rho
         y_coor = xi_rho
+    elif product=='NORKYST800':
+        x, y = find_nearest_cartCoord(ds.lon, ds.lat, lon, lat)
+        lon_near = ds.lon.sel(Y=y, X=x).values[0][0]
+        lat_near = ds.lat.sel(Y=y, X=x).values[0][0]  
+        x_coor = x
+        y_coor = y
     print('Found nearest: lon.='+str(lon_near)+',lat.=' + str(lat_near))     
     return x_coor, y_coor, lon_near, lat_near
 
@@ -176,6 +195,30 @@ def create_dataframe(product,ds, lon_near, lat_near,outfile,variable, start_time
         ds = ds.drop_vars('projection_lambert')
         ds = ds.drop_vars('latitude')
         ds = ds.drop_vars('longitude')
+    elif product=='NORKYST800': 
+        ds0 = ds
+        if 'depth' in ds['zeta'].dims:
+            ds['zeta'] = ds.zeta.sel(depth=0)
+
+        var_list = []
+        for var_name in variable:
+        # Check if 'depth' is not in the dimensions of the variable
+            if 'depth' in ds[var_name].dims:
+            # Append variable name to the list
+                var_list.append(var_name)
+
+        for i in range(len(height)):
+                variable_height = [k + '_'+str(height[i])+'m' for k in var_list]
+                ds[variable_height] = ds0[var_list].sel(depth=height[i],method='nearest')
+        ds = ds.drop_vars(var_list)
+        ds = ds.drop_vars('depth')
+        ds = ds.drop_vars('lat')
+        ds = ds.drop_vars('lon')
+        ds = ds.drop_vars('X')
+        ds = ds.drop_vars('Y')
+        ds = ds.drop_vars('projection_stere')
+        ds = ds.squeeze(drop=True)
+
     else: 
         drop_var = drop_variables(product=product) 
         ds = ds.drop_vars(drop_var)
@@ -184,7 +227,7 @@ def create_dataframe(product,ds, lon_near, lat_near,outfile,variable, start_time
     df = ds.to_dataframe()
     df = df.astype(float).round(2)
     df.index = pd.DatetimeIndex(data=ds.time.values)
-    
+
     top_header = '#'+product + ';LONGITUDE:'+str(lon_near.round(4))+';LATITUDE:' + str(lat_near.round(4)) 
     list_vars = [i for i in ds.data_vars]
     vars_info = ['#Variable_name;standard_name;long_name;units']
@@ -238,5 +281,25 @@ def read_commented_lines(datafile):
                 commented_lines = np.append(commented_lines,line)
     return commented_lines
 
+def remove_dimensions_from_netcdf(input_file, dimensions_to_remove=['X', 'Y']):
+    """
+    Remove specified dimensions from all variables in a NetCDF file and save the result to the same file.
+    
+    Parameters:
+    input_file (str): Path to the input NetCDF file.
+    dimensions_to_remove (list): A list of dimensions to remove from each variable.
+    """
+    # Check if the input file exists
+    if not os.path.exists(input_file):
+        return
 
-
+    # Open the input NetCDF file using xarray
+    ds = xr.open_dataset(input_file)
+
+    # Squeeze out the dimensions specified if they are singleton dimensions
+    for dim in dimensions_to_remove:
+        if dim in ds.dims and ds.sizes[dim] == 1:
+            ds = ds.squeeze(dim=dim)
+
+    # Write the updated dataset to the same NetCDF file
+    ds.to_netcdf(input_file)

metocean_api/ts/read_metno.py

@@ -8,8 +8,7 @@
 from nco import Nco
 from pathlib import Path
 
-
-from .aux_funcs import get_date_list, get_url_info, get_near_coord, create_dataframe, check_datafile_exists, read_commented_lines
+from .aux_funcs import *
 
 def NORAC_ts(self, save_csv = False, save_nc = False):
     """
@@ -204,6 +203,44 @@ def NORA3_stormsurge_ts(self, save_csv = False,save_nc = False):
 
     return df
 
+def NORKYST800_ts(self, save_csv = False, save_nc = False):
+    """
+    Extract times series of  the nearest gird point (lon,lat) from
+    Norkyst800.
+    """
+    self.variable.append('lon') # keep info of regular lon
+    self.variable.append('lat')  # keep info of regular lat
+    date_list = get_date_list(product=self.product, start_date=self.start_time, end_date=self.end_time)        
+
+    tempfile = tempfile_dir(self.product,self.lon, self.lat, date_list,dirName="cache")
+
+    # extract point and create temp files
+    for i in range(len(date_list)):
+        x_coor_str, y_coor_str, infile = get_url_info(product=self.product, date=date_list[i])
+
+        if i==0 or date_list[i].strftime('%Y-%m-%d %H:%M:%S') == '2019-02-27 00:00:00': # '2019-02-27' change to new model set up
+            x_coor, y_coor, lon_near, lat_near = get_near_coord(infile=infile, lon=self.lon, lat=self.lat, product=self.product)        
+
+        opt = ['-O -v '+",".join(self.variable)+' -d '+x_coor_str+','+str(x_coor.values[0])+' -d '+y_coor_str+','+str(y_coor.values[0])]
+
+        apply_nco(infile,tempfile[i],opt)
+        remove_dimensions_from_netcdf(tempfile[i], dimensions_to_remove=['X', 'Y'])
+
+    check_datafile_exists(self.datafile)
+    #merge temp files
+    existing_files = [f for f in tempfile if os.path.exists(f)]
+    ds = xr.open_mfdataset(paths=existing_files[:])
+    #Save in csv format   
+    df = create_dataframe(product=self.product,ds=ds, lon_near=lon_near, lat_near=lat_near, outfile=self.datafile, variable=self.variable[:-2], start_time = self.start_time, end_time = self.end_time, save_csv=save_csv,save_nc = save_nc, height=self.height)    
+    ds.close()
+    #remove temp files
+    #for i in range(len(date_list)):
+    #    os.remove(tempfile[i])
+
+    return df
+
+
+
 def NORA3_combined_ts(self, save_csv = True,save_nc = False):
     self.variable = ['hs','tp','fpI', 'tm1','tm2','tmp','Pdir','thq', 'hs_sea','tp_sea','thq_sea' ,'hs_swell','tp_swell','thq_swell']
     self.product = 'NORA3_wave_sub'

metocean_api/ts/ts_mod.py

@@ -93,6 +93,12 @@ def import_data(self, save_csv = True, save_nc = False):
       self.height = [50, 100, 150, 200, 300]
       self.variable =  ['wind_speed', 'wind_direction', 'air_temperature', 'relative_humidity', 'density', 'tke'] 
       self.data = NORA3_atm3hr_ts(self, save_csv = save_csv,save_nc = save_nc) 
+    elif self.product == 'NORKYST800':
+      self.height = [ 0.,    3.,   10.,   15.,   25.,   50.,   75.,  100.,  150.,  200.,
+        250.,  300.,  500., 1000., 2000., 3000.]
+      if self.variable == []:
+        self.variable =  ['salinity','temperature', 'u','v','zeta'] 
+      self.data = NORKYST800_ts(self, save_csv = save_csv,save_nc = save_nc) 
     elif self.product.startswith('E39'):
       if self.variable == []:
         self.variable =  ['Hm0'] 

tests/test_extract_data.py

@@ -59,3 +59,15 @@ def test_extract_OBS():
         pass
     else:
         raise ValueError("Shape is not correct")    
+
+def test_NORKYST800():
+    # Define TimeSeries-object
+    df_ts = ts.TimeSeries(lon=3.73, lat=64.60,start_time='2020-09-14', end_time='2020-09-14', product='NORKYST800')
+    # Import data from thredds.met.no 
+    df_ts.import_data(save_csv=False,save_nc=False)
+    if df_ts.data.shape == (24, 65):
+        pass
+    else:
+        raise ValueError("Shape is not correct")   
+
+

version.py

@@ -1,4 +1,4 @@
-__version__ = "1.1.2"
+__version__ = "1.1.3"
 
 
 def git_describe():