Initial commit

00_download_ERA5_ivt.py

@@ -0,0 +1,86 @@
+# How to use the CDS API:
+# https://cds.climate.copernicus.eu/api-how-to
+
+import os
+import argparse
+import yaml
+
+import cdsapi
+
+# argument parameters
+parser = argparse.ArgumentParser(
+    description=__doc__,
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+)
+parser.add_argument(
+    '--config', '-c',
+    type=str,
+    help='path to the config.yml file',
+    default=os.path.join(os.getcwd(), 'config.yml'),
+)
+args = parser.parse_args()
+
+# load config
+config = yaml.safe_load(open(args.config))
+
+# parameters
+years = range(config['year_start'], config['year_end'] + 1)
+variables = [
+    'vertical_integral_of_eastward_water_vapour_flux',
+    'vertical_integral_of_northward_water_vapour_flux',
+]
+
+# filesystem
+data_folder_e = config['data_folder_e']
+data_folder_n = config['data_folder_n']
+
+# initiate client
+c = cdsapi.Client()
+
+# get data
+for year in years:
+    for variable in variables:
+
+        if "eastward" in variable:
+            folder = data_folder_e
+        elif "northward" in variable:
+            folder = data_folder_n
+
+        c.retrieve(
+            'reanalysis-era5-single-levels',
+            {
+                'product_type': 'reanalysis',
+                'format': 'netcdf',
+                'variable': variable,
+                'year': year,
+                'month': [
+                    '01', '02', '03',
+                    '04', '05', '06',
+                    '07', '08', '09',
+                    '10', '11', '12',
+                ],
+                'day': [
+                    '01', '02', '03',
+                    '04', '05', '06',
+                    '07', '08', '09',
+                    '10', '11', '12',
+                    '13', '14', '15',
+                    '16', '17', '18',
+                    '19', '20', '21',
+                    '22', '23', '24',
+                    '25', '26', '27',
+                    '28', '29', '30',
+                    '31',
+                ],
+                'time': [
+                    '00:00', '01:00', '02:00',
+                    '03:00', '04:00', '05:00',
+                    '06:00', '07:00', '08:00',
+                    '09:00', '10:00', '11:00',
+                    '12:00', '13:00', '14:00',
+                    '15:00', '16:00', '17:00',
+                    '18:00', '19:00', '20:00',
+                    '21:00', '22:00', '23:00',
+                ],
+            },
+            os.path.join(folder, f'{year}.nc'))

01_regrid_ivt.py

@@ -0,0 +1,132 @@
+# Copyright (C) 2022 by
+# Dominik Traxl <dominik.traxl@posteo.org>
+# All rights reserved.
+# GPL-3.0 license.
+
+import os
+import argparse
+
+import yaml
+import numpy as np
+import pandas as pd
+import xarray as xr
+import xesmf as xe
+
+# argument parameters
+parser = argparse.ArgumentParser(
+    description=__doc__,
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+)
+parser.add_argument(
+    'year',
+    type=str,
+    help='regrid given year',
+)
+parser.add_argument(
+    '--config', '-c',
+    type=str,
+    help='path to the config.yml file',
+    default=os.path.join(os.getcwd(), 'config.yml'),
+)
+args = parser.parse_args()
+
+# load config
+config = yaml.safe_load(open(args.config))
+
+# config parameters
+spatial_resolution = config['spatial_resolution']
+temporal_resolution = config['temporal_resolution']
+
+# filesystem
+data_folder_e = config['data_folder_e']
+data_folder_n = config['data_folder_n']
+output_folder = config['output_folder']
+os.makedirs(
+    os.path.join(output_folder, 'ivt_regridded', args.year),
+    exist_ok=True
+)
+
+# era5 files
+# name: vertical integral of eastward water vapour flux; yearly downloads
+vapfut_files = os.listdir(data_folder_e)
+# name: vertical integral of northward water vapour flux; yearly downloads
+vapfvt_files = os.listdir(data_folder_n)
+
+# output resolution
+ds_out = xr.Dataset({
+    "latitude": (["latitude"],
+                 np.arange(90, -90-spatial_resolution, -spatial_resolution)),
+    "longitude": (["longitude"],
+                  np.arange(0, 360, spatial_resolution))
+})
+
+# create file dictionary
+uflux_files = {}
+for f in vapfut_files:
+    year = f.split('_')[-1][:4]
+    uflux_files[year] = f
+
+vflux_files = {}
+for f in vapfvt_files:
+    year = f.split('_')[-1][:4]
+    vflux_files[year] = f
+
+# load uflux
+uflux_fn = uflux_files[args.year]
+ds_uflux = xr.open_dataset(os.path.join(data_folder_e, uflux_fn))
+assert len(ds_uflux.data_vars) == 1
+ds_uflux = ds_uflux.rename({list(ds_uflux.data_vars)[0]: 'uflux'})
+
+# load vflux
+vflux_fn = vflux_files[args.year]
+ds_vflux = xr.open_dataset(os.path.join(data_folder_n, vflux_fn))
+assert len(ds_vflux.data_vars) == 1
+ds_vflux = ds_vflux.rename({list(ds_vflux.data_vars)[0]: 'vflux'})
+
+# combine
+ds = ds_uflux
+ds['vflux'] = ds_vflux['vflux']
+
+# assert correct number of time-slices
+days_in_year = pd.Timestamp(int(args.year), 12, 31).day_of_year
+assert ds.time.shape[0] == days_in_year * 24
+
+# pos array
+positions = np.arange(0, ds.time.shape[0] + 24, 24)
+
+# regrid on a daily basis
+for i in range(len(positions) - 1):
+
+    print(f'{args.year}-{i+1:03d}')
+
+    # subset
+    dst = ds.isel(time=slice(positions[i], positions[i+1]))
+
+    # compute ivt
+    dst['ivt'] = np.sqrt(dst['uflux']**2 + dst['vflux']**2)
+
+    # create regridder
+    try:
+        regridder = xe.Regridder(
+            dst, ds_out, "bilinear", periodic=True,
+            reuse_weights=True,
+            filename=os.path.join(output_folder,
+                                  "ivt_regridded",
+                                  "bilinear_peri.nc"))
+    except OSError as e:
+        print(OSError, e)
+        regridder = xe.Regridder(dst, ds_out, "bilinear", periodic=True)
+        regridder.to_netcdf(
+            os.path.join(output_folder, 'ivt_regridded', 'bilinear_peri.nc'))
+
+    # perform regridding
+    dst_r = regridder(dst)
+
+    # resample time
+    dst_rr = dst_r.resample(time=f'{temporal_resolution}h',
+                            loffset='3h',
+                            ).mean()
+
+    # store
+    dst_rr.to_netcdf(os.path.join(
+        output_folder, 'ivt_regridded', args.year, f'{i:03d}.nc'), mode='w')

02_aggregate.py

@@ -0,0 +1,42 @@
+# Copyright (C) 2022 by
+# Dominik Traxl <dominik.traxl@posteo.org>
+# All rights reserved.
+# GPL-3.0 license.
+
+import os
+import sys
+import argparse
+
+import yaml
+import xarray as xr
+
+# argument parameters
+parser = argparse.ArgumentParser(
+    description=__doc__,
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+)
+parser.add_argument(
+    '--config', '-c',
+    type=str,
+    help='path to the config.yml file',
+    default=os.path.join(os.getcwd(), 'config.yml'),
+)
+args = parser.parse_args()
+
+# load config
+config = yaml.safe_load(open(args.config))
+
+# exit if user does not want to regrid
+if not config['do_regridding']:
+    sys.exit(0)
+
+# filesystem
+output_folder = config['output_folder']
+
+# aggregate files for each year
+for year in range(config['year_start'], config['year_end'] + 1):
+    print(f'opening {year}')
+    ds = xr.open_mfdataset(os.path.join(
+        output_folder, 'ivt_regridded', str(year), '*.nc'))
+    print('writing to file')
+    ds.to_netcdf(os.path.join(output_folder, 'ivt_regridded', f'{year}.nc'))

03_ipart_ar_tracking_thr_multifile.py

@@ -0,0 +1,55 @@
+# Copyright (C) 2022 by
+# Dominik Traxl <dominik.traxl@posteo.org>
+# All rights reserved.
+# GPL-3.0 license.
+
+import os
+import argparse
+
+import yaml
+from ipart import thr
+
+# argument parameters
+parser = argparse.ArgumentParser(
+    description=__doc__,
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+)
+parser.add_argument(
+    '--config', '-c',
+    type=str,
+    help='path to the config.yml file',
+    default=os.path.join(os.getcwd(), 'config.yml'),
+)
+args = parser.parse_args()
+
+# load config
+config = yaml.safe_load(open(args.config))
+
+# parameters
+years = range(config['year_start'], config['year_end'] + 1)
+kernel = config['kernel']
+shift_lon = config['shift_lon']
+
+# filesystem
+output_folder = config['output_folder']
+os.makedirs(os.path.join(output_folder, 'ipart', 'thr'), exist_ok=True)
+if config['do_regridding']:
+    input_ivt_folder = 'ivt_regridded/'
+else:
+    input_ivt_folder = 'ivt/'
+
+# create file list
+filelist = []
+for year in years:
+    file_in_name = f'{year:d}.nc'
+    abpath_in = os.path.join(output_folder, input_ivt_folder, file_in_name)
+    filelist.append(abpath_in)
+print(filelist)
+
+# need at least two years of data
+assert len(filelist) >= 2
+
+# compute
+thr.rotatingTHR(filelist, 'ivt', kernel,
+                os.path.join(output_folder, 'ipart', 'thr'),
+                shift_lon=shift_lon)