From 598bb4da557bc850b58e72db801e14c281d8034c Mon Sep 17 00:00:00 2001 From: KonstantinChri Date: Sun, 23 Jun 2024 01:36:10 +0200 Subject: [PATCH 1/4] add norkyst800 --- docs/index.rst | 3 +++ metocean_api/ts/aux_funcs.py | 34 +++++++++++++++++++++++++++++++-- metocean_api/ts/read_metno.py | 36 +++++++++++++++++++++++++++++++++++ metocean_api/ts/ts_mod.py | 6 ++++++ tests/test_extract_data.py | 10 ++++++++++ version.py | 2 +- 6 files changed, 88 insertions(+), 3 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index 653c4ad..dd2dc20 100755 --- a/docs/index.rst +++ b/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 (temperature, currents, salinity etc) Norkyst800 data: product='NORKYST800' + + Dataset: https://thredds.met.no/thredds/fou-hi/norkyst800v2.html * For global reananalysis ERA5 (wind and waves): product='ERA5' diff --git a/metocean_api/ts/aux_funcs.py b/metocean_api/ts/aux_funcs.py index 8546ea1..d19d14f 100644 --- a/metocean_api/ts/aux_funcs.py +++ b/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,11 @@ 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': + infile = 'https://thredds.met.no/thredds/dodsC/fou-hi/norkyst800m-1h/NorKyst-800m_ZDEPTHS_his.an.'+date.strftime('%Y%m%d%H')+'.nc' + #https://thredds.met.no/thredds/dodsC/fou-hi/norkyst800m-1h/NorKyst-800m_ZDEPTHS_his.an.2024062200.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 +99,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 +122,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 +155,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 +193,19 @@ 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 + for i in range(len(height)): + variable_height = [k + '_'+str(height[i])+'m' for k in variable] + ds[variable_height] = ds0[variable].sel(depth=height[i]) + ds = ds.drop_vars(variable) + 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) diff --git a/metocean_api/ts/read_metno.py b/metocean_api/ts/read_metno.py index 61843de..0f6a9a2 100644 --- a/metocean_api/ts/read_metno.py +++ b/metocean_api/ts/read_metno.py @@ -204,6 +204,42 @@ 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: + 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) + + check_datafile_exists(self.datafile) + #merge temp files + ds = xr.open_mfdataset(paths=tempfile[:]) + #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' diff --git a/metocean_api/ts/ts_mod.py b/metocean_api/ts/ts_mod.py index f2254aa..fea7e0c 100644 --- a/metocean_api/ts/ts_mod.py +++ b/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','tke','temperature', 'u','v','w','zeta', 'ubar', 'vbar'] + self.data = NORKYST800_ts(self, save_csv = save_csv,save_nc = save_nc) elif self.product.startswith('E39'): if self.variable == []: self.variable = ['Hm0'] diff --git a/tests/test_extract_data.py b/tests/test_extract_data.py index 0c44f29..4d4959b 100644 --- a/tests/test_extract_data.py +++ b/tests/test_extract_data.py @@ -59,3 +59,13 @@ def test_extract_OBS(): pass else: raise ValueError("Shape is not correct") + +def test_NORKYST800(): + # Define TimeSeries-object + df_ts = ts.TimeSeries(lon=1.320, lat=53.324,start_time='2020-01-01', end_time='2020-01-02', product='NORKYST800') + # Import data from thredds.met.no + df_ts.import_data(save_csv=False,save_nc=False) + if df_ts.data.shape == (48, 144): + pass + else: + raise ValueError("Shape is not correct") diff --git a/version.py b/version.py index f6f23ec..c2fe38a 100644 --- a/version.py +++ b/version.py @@ -1,4 +1,4 @@ -__version__ = "1.1.2" +__version__ = "1.1.3" def git_describe(): From 9a33c9c2950b2fc8efe1f3eed9d8627d4573c728 Mon Sep 17 00:00:00 2001 From: KonstantinChri Date: Tue, 25 Jun 2024 14:32:25 +0200 Subject: [PATCH 2/4] update --- metocean_api/ts/aux_funcs.py | 12 ++++++++---- metocean_api/ts/ts_mod.py | 2 +- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/metocean_api/ts/aux_funcs.py b/metocean_api/ts/aux_funcs.py index d19d14f..9288de6 100644 --- a/metocean_api/ts/aux_funcs.py +++ b/metocean_api/ts/aux_funcs.py @@ -77,8 +77,10 @@ def get_url_info(product, date): x_coor_str = 'eta_rho' y_coor_str = 'xi_rho' elif product == 'NORKYST800': - infile = 'https://thredds.met.no/thredds/dodsC/fou-hi/norkyst800m-1h/NorKyst-800m_ZDEPTHS_his.an.'+date.strftime('%Y%m%d%H')+'.nc' - #https://thredds.met.no/thredds/dodsC/fou-hi/norkyst800m-1h/NorKyst-800m_ZDEPTHS_his.an.2024062200.nc + if date Date: Sat, 13 Jul 2024 13:30:33 +0200 Subject: [PATCH 3/4] fix NORKYST800 old and new model set up --- docs/index.rst | 2 +- examples/example_import_NORA3.py | 2 +- metocean_api/ts/aux_funcs.py | 46 +++++++++++++++++++++++++------- metocean_api/ts/read_metno.py | 11 ++++---- tests/test_extract_data.py | 6 +++-- 5 files changed, 48 insertions(+), 19 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index dd2dc20..b91a940 100755 --- a/docs/index.rst +++ b/docs/index.rst @@ -96,7 +96,7 @@ 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 (temperature, currents, salinity etc) Norkyst800 data: product='NORKYST800' +* 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 diff --git a/examples/example_import_NORA3.py b/examples/example_import_NORA3.py index 34a7af6..2666534 100644 --- a/examples/example_import_NORA3.py +++ b/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) diff --git a/metocean_api/ts/aux_funcs.py b/metocean_api/ts/aux_funcs.py index 9288de6..55ff721 100644 --- a/metocean_api/ts/aux_funcs.py +++ b/metocean_api/ts/aux_funcs.py @@ -77,9 +77,9 @@ def get_url_info(product, date): 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' - else: + 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' @@ -197,12 +197,20 @@ def create_dataframe(product,ds, lon_near, lat_near,outfile,variable, start_time ds = ds.drop_vars('longitude') elif product=='NORKYST800': ds0 = ds - breakpoint() + 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)): - if 'depth' in ds0[variable].dims: # Check if 'depth' is a dimension of the variable - variable_height = [k + '_'+str(height[i])+'m' for k in variable] - ds[variable_height] = ds0[variable].sel(depth=height[i]) - ds = ds.drop_vars(variable) + 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') @@ -210,6 +218,7 @@ def create_dataframe(product,ds, lon_near, lat_near,outfile,variable, start_time 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) @@ -218,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'] @@ -272,5 +281,22 @@ 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 a new file. + + Parameters: + input_file (str): Path to the input NetCDF file. + output_file (str): Path to the output NetCDF file. + dimensions_to_remove (list): A list of dimensions to remove from each variable. + """ + # 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 output NetCDF file + ds.to_netcdf(input_file) diff --git a/metocean_api/ts/read_metno.py b/metocean_api/ts/read_metno.py index 0f6a9a2..4a6cc91 100644 --- a/metocean_api/ts/read_metno.py +++ b/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): """ @@ -218,17 +217,19 @@ def NORKYST800_ts(self, save_csv = False, save_nc = False): # 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: + + 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 - ds = xr.open_mfdataset(paths=tempfile[:]) + 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() diff --git a/tests/test_extract_data.py b/tests/test_extract_data.py index 4d4959b..f50a89c 100644 --- a/tests/test_extract_data.py +++ b/tests/test_extract_data.py @@ -62,10 +62,12 @@ def test_extract_OBS(): def test_NORKYST800(): # Define TimeSeries-object - df_ts = ts.TimeSeries(lon=1.320, lat=53.324,start_time='2020-01-01', end_time='2020-01-02', product='NORKYST800') + 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 == (48, 144): + if df_ts.data.shape == (24, 65): pass else: raise ValueError("Shape is not correct") + + From 6fff993df292262a4efca4605bee567d65b0d5c1 Mon Sep 17 00:00:00 2001 From: KonstantinChri Date: Sat, 13 Jul 2024 14:28:17 +0200 Subject: [PATCH 4/4] Check if the input file exists in remove_dimensions_from_netcdf --- metocean_api/ts/aux_funcs.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/metocean_api/ts/aux_funcs.py b/metocean_api/ts/aux_funcs.py index 55ff721..02c1056 100644 --- a/metocean_api/ts/aux_funcs.py +++ b/metocean_api/ts/aux_funcs.py @@ -283,13 +283,16 @@ def read_commented_lines(datafile): 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 a new file. + 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. - output_file (str): Path to the output 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) @@ -298,5 +301,5 @@ def remove_dimensions_from_netcdf(input_file, dimensions_to_remove=['X', 'Y']): if dim in ds.dims and ds.sizes[dim] == 1: ds = ds.squeeze(dim=dim) - # Write the updated dataset to the output NetCDF file - ds.to_netcdf(input_file) + # Write the updated dataset to the same NetCDF file + ds.to_netcdf(input_file) \ No newline at end of file