From 09c45e7d91e418283a5e70c67c183ebba1a7b1c7 Mon Sep 17 00:00:00 2001 From: NabeelAI Date: Sat, 4 Apr 2020 15:45:24 +0500 Subject: [PATCH 1/8] fix(): discrepencies in netcdf export --- pysteps/io/exporters.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/pysteps/io/exporters.py b/pysteps/io/exporters.py index f70bf15..699effd 100644 --- a/pysteps/io/exporters.py +++ b/pysteps/io/exporters.py @@ -413,20 +413,16 @@ def initialize_forecast_exporter_netcdf(filename, startdate, timestep, var_time = ncf.createVariable("time", np.int, dimensions=("time",)) if incremental != "timestep": - if product == 'precip_probability': - var_time[:] = [i*timestep for i in range(1, n_timesteps+1)] - else: - var_time[:] = [i*timestep*60 for i in range(1, n_timesteps+1)] + var_time[:] = [i*timestep for i in range(1, n_timesteps+1)] var_time.long_name = "forecast time" startdate_str = datetime.strftime(startdate, "%Y-%m-%d %H:%M:%S") - var_time.units = "minutes since %s" % startdate_str if product == 'precip_probability' \ - else "seconds since %s" % startdate_str - - dimensions = ("time", "y", "x") + # var_time.units = "minutes since %s" % startdate_str if product == 'precip_probability' \ + # else "seconds since %s" % startdate_str + var_time.units = "minutes since %s" % startdate_str var_F = ncf.createVariable(var_name, np.float32, - dimensions=dimensions, + dimensions=("time", "y", "x"), zlib=True, complevel=9) if var_standard_name is not None: From 59c5f94e58f5c2c2254362377a05cc2d136c23a4 Mon Sep 17 00:00:00 2001 From: NabeelAI Date: Sun, 5 Apr 2020 18:08:16 +0500 Subject: [PATCH 2/8] fix(): fix netcdf export issue --- examples/probabilistic_nowcast_with_export.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/probabilistic_nowcast_with_export.py b/examples/probabilistic_nowcast_with_export.py index 47b52a4..f306bad 100755 --- a/examples/probabilistic_nowcast_with_export.py +++ b/examples/probabilistic_nowcast_with_export.py @@ -177,7 +177,7 @@ # set -1 for nan prob_array[np.isnan(prob_array)] = -1 export_initializer = stp.io.get_method('netcdf', 'exporter') -exporter = export_initializer(filename, startdate, timestep, n_leadtimes , shape, n_ens_members, metadata, +exporter = export_initializer(filename, startdate, timestep, n_leadtimes, shape, n_ens_members, metadata, product='precip_probability', incremental=None) stp.io.export_forecast_dataset(prob_array, exporter) stp.io.close_forecast_file(exporter) From f7d6064834923688ac23083b3af1e8976a258d25 Mon Sep 17 00:00:00 2001 From: NabeelAI Date: Sun, 19 Apr 2020 21:43:35 +0500 Subject: [PATCH 3/8] feat(): allow one minute forecast in precipitation nowcasting --- examples/precipitation_nowcast.py | 53 ++++++++++++++++++++----------- pysteps/nowcasts/extrapolation.py | 3 +- 2 files changed, 36 insertions(+), 20 deletions(-) diff --git a/examples/precipitation_nowcast.py b/examples/precipitation_nowcast.py index 3f003dd..6ca5427 100644 --- a/examples/precipitation_nowcast.py +++ b/examples/precipitation_nowcast.py @@ -2,6 +2,9 @@ # -*- coding: utf-8 -*- import os, sys import pyximport + +from pysteps_custom_utils.utils import convert_dbz_to_mm + pyximport.install() from tendo import singleton @@ -68,7 +71,7 @@ ## forecast parameters n_prvs_times = 3 # use at least 9 with DARTS -n_leadtimes = 5 +n_leadtimes = 40 n_ens_members = 3 n_cascade_levels = 8 ar_order = 2 @@ -82,6 +85,7 @@ transformation = "dB" # None or dB adjust_domain = None # None or square seed = 42 # for reproducibility +output_time_step = 1 vel_pert_kwargs = dict() vel_pert_kwargs["p_par"] = [ 0.38550792, 0.62097167, -0.23937287] @@ -91,7 +95,7 @@ # Read-in the data print('Read the data...', startdate_str) startdate = datetime.datetime.strptime(startdate_str, "%Y%m%d%H%M") - +# startdate = datetime.datetime(2019, 10, 7, 18, 40) ## import data specifications ds = cfg.get_specifications(data_source) @@ -103,6 +107,9 @@ importer = stp.io.get_method(ds.importer, "importer") R, _, metadata = stp.io.read_timeseries(input_files, importer, **ds.importer_kwargs) +# Preprocess the data +# R = convert_dbz_to_mm(R) + # Prepare input files print("Prepare the data...") @@ -133,6 +140,10 @@ oflow_method = stp.motion.get_method(oflow_method) UV = oflow_method(R) +# Change motion field according to timeStep +if output_time_step < 10: + UV = UV - (UV/10)*(10-output_time_step) + # Perform the nowcast extrap_kwargs={} extrap_kwargs['allow_nonfinite_values'] = True @@ -140,22 +151,26 @@ # apply nan mask back R[nan_mask] = np.nan -nwc_method = stp.nowcasts.get_method(nwc_method) -R_fct = nwc_method(R[-3:, :, :], UV, - n_leadtimes, - n_ens_members, - n_cascade_levels=8, - R_thr=-10.0, - kmperpixel=metadata["xpixelsize"]/1000, - timestep=5, - decomp_method="fft", - bandpass_filter_method="gaussian", - noise_method="nonparametric", - vel_pert_method="bps", - mask_method="incremental", - seed=seed, - conserve_radar_mask=True, - vel_pert_kwargs=vel_pert_kwargs)[-1, :, :] +extrapolate = stp.nowcasts.get_method("extrapolation") +R_fct = extrapolate(R[-1], UV, n_leadtimes, allow_nans=True) + + +# nwc_method = stp.nowcasts.get_method(nwc_method) +# R_fct = nwc_method(R[-3:, :, :], UV, +# n_leadtimes, +# n_ens_members, +# n_cascade_levels=8, +# R_thr=-10.0, +# kmperpixel=metadata["xpixelsize"]/1000, +# timestep=output_time_step, +# decomp_method="fft", +# bandpass_filter_method="gaussian", +# noise_method="nonparametric", +# vel_pert_method="bps", +# mask_method="incremental", +# seed=seed, +# conserve_radar_mask=True, +# vel_pert_kwargs=vel_pert_kwargs)[-1, :, :] ## convert all data to mm/h converter = stp.utils.get_method("mm/h") @@ -177,7 +192,7 @@ # # set -1 for nan R_fct[np.isnan(R_fct)] = -1 export_initializer = stp.io.get_method('netcdf', 'exporter') -exporter = export_initializer(filename, startdate, timestep, n_leadtimes, shape, 1, metadata, +exporter = export_initializer(filename, startdate, output_time_step, n_leadtimes, shape, 1, metadata, incremental=None) stp.io.export_forecast_dataset(R_fct, exporter) stp.io.close_forecast_file(exporter) diff --git a/pysteps/nowcasts/extrapolation.py b/pysteps/nowcasts/extrapolation.py index f0dff8a..bc6971a 100644 --- a/pysteps/nowcasts/extrapolation.py +++ b/pysteps/nowcasts/extrapolation.py @@ -17,7 +17,7 @@ def forecast(precip, velocity, num_timesteps, extrap_method="semilagrangian", extrap_kwargs=None, - measure_time=False): + measure_time=False, allow_nans=False): """Generate a nowcast by applying a simple advection-based extrapolation to the given precipitation field. @@ -73,6 +73,7 @@ def forecast(precip, velocity, num_timesteps, extrapolation_method = extrapolation.get_method(extrap_method) precip_forecast = extrapolation_method(precip, velocity, num_timesteps, + allow_nonfinite_values=allow_nans, **extrap_kwargs) if measure_time: From 850ca158be8d8b6eefbd15ae58612ba3fa58fdd0 Mon Sep 17 00:00:00 2001 From: Nabeel07925 Date: Sun, 3 May 2020 17:47:16 +0500 Subject: [PATCH 4/8] feat(): add script for precipitation accumulation export with one hour sum --- .../precipitation_accumulation_nowcast.py | 204 ++++++++++++++++++ 1 file changed, 204 insertions(+) diff --git a/examples/precipitation_accumulation_nowcast.py b/examples/precipitation_accumulation_nowcast.py index e69de29..a3bc6e5 100644 --- a/examples/precipitation_accumulation_nowcast.py +++ b/examples/precipitation_accumulation_nowcast.py @@ -0,0 +1,204 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +import os, sys +import pyximport +import copy + +from pysteps_custom_utils.utils import convert_dbz_to_mm + +pyximport.install() + +from tendo import singleton + +me = singleton.SingleInstance() # will sys.exit(-1) if other instance is running + + +"""Stochastic ensemble precipitation nowcasting + +The script shows how to run a stochastic ensemble of precipitation nowcasts with +pysteps. + +More info: https://pysteps.github.io/ +""" +import datetime +import numpy as np +sys.path.append(os.path.join(os.path.dirname(__file__), '../')) + +import pysteps as stp +import config as cfg + +# List of case studies that can be used in this tutorial + +#+-------+--------------+-------------+----------------------------------------+ +#| event | start_time | data_source | description | +#+=======+==============+=============+========================================+ +#| 01 | 201701311030 | mch | orographic precipitation | +#+-------+--------------+-------------+----------------------------------------+ +#| 02 | 201505151630 | mch | non-stationary field, apparent rotation| +#+-------+--------------+------------------------------------------------------+ +#| 03 | 201609281530 | fmi | stratiform rain band | +#+-------+--------------+-------------+----------------------------------------+ +#| 04 | 201705091130 | fmi | widespread convective activity | +#+-------+--------------+-------------+----------------------------------------+ +#| 05 | 201806161100 | bom | bom example data | +#+-------+--------------+-------------+----------------------------------------+ + +# Set parameters for this tutorial +data_source="gimet" + +## import data specifications +ds = cfg.get_specifications(data_source) + + +## input data (copy/paste values from table above) +archive_dir=ds.root_path+"/"+max(os.listdir(ds.root_path)) +last_fname=max(os.listdir(archive_dir)) +print(last_fname) +startdate_str=last_fname[-18:-10]+last_fname[-9:-5] + +print(startdate_str) + +## methods +oflow_method = "lucaskanade" # lucaskanade, darts, None +nwc_method = "steps" +adv_method = "semilagrangian" # semilagrangian, eulerian +noise_method = "nonparametric" # parametric, nonparametric, ssft +bandpass_filter = "gaussian" +decomp_method = "fft" + +## forecast parameters +n_prvs_times = 6 # use at least 9 with DARTS +n_leadtimes = 6 +r_threshold = 0.1 # rain/no-rain threshold [mm/h] +unit = "mm/h" # mm/h or dBZ +transformation = "dB" # None or dB +adjust_domain = None # None or square +output_time_step = 1 + +vel_pert_kwargs = dict() +vel_pert_kwargs["p_par"] = [ 0.38550792, 0.62097167, -0.23937287] +vel_pert_kwargs["p_perp"] = [0.2240485, 0.68900218, 0.24242502] + + +# Read-in the data +print('Read the data...', startdate_str) +startdate = datetime.datetime.strptime(startdate_str, "%Y%m%d%H%M") +# startdate = datetime.datetime(2019, 10, 7, 18, 40) +## import data specifications +ds = cfg.get_specifications(data_source) + +## find radar field filenames +input_files = stp.io.find_by_date(startdate, ds.root_path, ds.path_fmt, ds.fn_pattern, + ds.fn_ext, ds.timestep, n_prvs_times, 0) + +## read radar field files +importer = stp.io.get_method(ds.importer, "importer") +R, _, metadata = stp.io.read_timeseries(input_files, importer, **ds.importer_kwargs) + +# Preprocess the data +# R = convert_dbz_to_mm(R) + +# Prepare input files +print("Prepare the data...") + +## if requested, make sure we work with a square domain +reshaper = stp.utils.get_method(adjust_domain) +R, metadata = reshaper(R, metadata, method="pad") + +## if necessary, convert to rain rates [mm/h] +converter = stp.utils.get_method("mm/h") +R, metadata = converter(R, metadata) + +## threshold the data +R[R Date: Sun, 3 May 2020 18:06:45 +0500 Subject: [PATCH 5/8] feat(): support 2 hours precipitation accumulation as well. --- examples/precipitation_accumulation_nowcast.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/precipitation_accumulation_nowcast.py b/examples/precipitation_accumulation_nowcast.py index a3bc6e5..3093c7a 100644 --- a/examples/precipitation_accumulation_nowcast.py +++ b/examples/precipitation_accumulation_nowcast.py @@ -68,7 +68,7 @@ ## forecast parameters n_prvs_times = 6 # use at least 9 with DARTS -n_leadtimes = 6 +n_leadtimes = 6 # use 6 for one hour extrapolation and 12 for 2 hours extrapolation r_threshold = 0.1 # rain/no-rain threshold [mm/h] unit = "mm/h" # mm/h or dBZ transformation = "dB" # None or dB @@ -166,7 +166,7 @@ print(R_computed.shape) # compute sum at one hour -R_final = np.zeros((R.shape[0], R.shape[1], R.shape[2])) +R_final = np.zeros((R_all.shape[0]-6, R.shape[1], R.shape[2])) for i in range(R_final.shape[0]): if i == 0: R_final[i] = np.sum(R_computed[:6], axis=0) From a83b9989c6fb0f22b990c0bb75b3cbc505f51bce Mon Sep 17 00:00:00 2001 From: Eugene Savelov Date: Mon, 11 May 2020 18:12:29 +0300 Subject: [PATCH 6/8] fix compile errors --- examples/precipitation_accumulation_nowcast.py | 2 +- pysteps/io/importers.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/precipitation_accumulation_nowcast.py b/examples/precipitation_accumulation_nowcast.py index 3093c7a..9daabf4 100644 --- a/examples/precipitation_accumulation_nowcast.py +++ b/examples/precipitation_accumulation_nowcast.py @@ -4,7 +4,7 @@ import pyximport import copy -from pysteps_custom_utils.utils import convert_dbz_to_mm +#from pysteps_custom_utils.utils import convert_dbz_to_mm pyximport.install() diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index 410ff6f..a894922 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -115,7 +115,7 @@ except ImportError: netcdf4_imported = False try: - import PIL + import PIL.Image pil_imported = True except ImportError: From 3c700abb255d7a5a6ae8bab4200db34cd602b713 Mon Sep 17 00:00:00 2001 From: Nabeel07925 Date: Sun, 17 May 2020 16:27:12 +0500 Subject: [PATCH 7/8] feat(): make a and b coefficients configurable --- examples/precipitation_accumulation_nowcast.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/examples/precipitation_accumulation_nowcast.py b/examples/precipitation_accumulation_nowcast.py index 3093c7a..35f1bc8 100644 --- a/examples/precipitation_accumulation_nowcast.py +++ b/examples/precipitation_accumulation_nowcast.py @@ -60,11 +60,6 @@ ## methods oflow_method = "lucaskanade" # lucaskanade, darts, None -nwc_method = "steps" -adv_method = "semilagrangian" # semilagrangian, eulerian -noise_method = "nonparametric" # parametric, nonparametric, ssft -bandpass_filter = "gaussian" -decomp_method = "fft" ## forecast parameters n_prvs_times = 6 # use at least 9 with DARTS @@ -74,11 +69,8 @@ transformation = "dB" # None or dB adjust_domain = None # None or square output_time_step = 1 - -vel_pert_kwargs = dict() -vel_pert_kwargs["p_par"] = [ 0.38550792, 0.62097167, -0.23937287] -vel_pert_kwargs["p_perp"] = [0.2240485, 0.68900218, 0.24242502] - +a = 316.0 +b = 1.5 # Read-in the data print('Read the data...', startdate_str) @@ -107,7 +99,7 @@ ## if necessary, convert to rain rates [mm/h] converter = stp.utils.get_method("mm/h") -R, metadata = converter(R, metadata) +R, metadata = converter(R, metadata, a=a, b=b) ## threshold the data R[R Date: Mon, 18 May 2020 00:22:02 +0500 Subject: [PATCH 8/8] fix(): fix accumulation --- examples/precipitation_accumulation_nowcast.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/precipitation_accumulation_nowcast.py b/examples/precipitation_accumulation_nowcast.py index eb80a01..fd78178 100644 --- a/examples/precipitation_accumulation_nowcast.py +++ b/examples/precipitation_accumulation_nowcast.py @@ -110,8 +110,8 @@ # R, metadata = converter(R, metadata) ## transform the data -transformer = stp.utils.get_method(transformation) -R, metadata = transformer(R, metadata) +# transformer = stp.utils.get_method(transformation) +# R, metadata = transformer(R, metadata) nan_mask = np.ma.masked_invalid(R).mask R[nan_mask] = metadata["zerovalue"] # to compute optical flow @@ -139,8 +139,8 @@ # R_all = np.ma.masked_less_equal(R_all, 0) # per minute precipitation -R_all[R_all<0] = 0 -R_all = R_all/60 +# R_all[R_all<0] = 0 +# R_all = R_all/60 slow_UV = UV - (UV/10)*(10-output_time_step) # one minute motion vectors R_computed = np.zeros((R_all.shape[0], R.shape[1], R.shape[2])) @@ -150,14 +150,16 @@ # compute sum at 10 min timestep for i in range(R_all.shape[0]): if i == 0: - R_computed[i] = R_all[i] + R_computed[i] = R_all[i]/60 ten_min_sum = R_computed[i] else: + prev_min_nowcast = R_all[i] for j in range(10): - one_min_nowcast = extrapolate(ten_min_sum, slow_UV, 1, allow_nans=True) - ten_min_sum += one_min_nowcast[0] + one_min_nowcast = extrapolate(prev_min_nowcast, slow_UV, 1, allow_nans=True) + ten_min_sum += one_min_nowcast[0]/60 + prev_min_nowcast = one_min_nowcast[0] R_computed[i] = ten_min_sum - ten_min_sum = R_all[i] + ten_min_sum = R_all[i]/60 print(R_computed.shape)