See also final script in ../scripts/freezing_level.py
In [1]:
%matplotlib inline
In [2]:
import sys
import os
aps_path = os.path.dirname(os.path.abspath("."))
if aps_path not in sys.path:
sys.path.append(aps_path)
print(aps_path, sys.path)
In [3]:
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
sns.set(style="dark")
import aps_io.get_arome as ga
from load_region import load_region, clip_region
In [4]:
hour_range = [24, 48]
ncfile = r"\\hdata\grid\metdata\prognosis\meps\det\archive\2019\meps_det_extracted_1km_20190307T00Z.nc"
#jd, altitude, land_area_fraction, nc_vars = ga.nc_load(ncfile, ["altitude_of_0_degree_isotherm"], time_period=[6, 25])
times, altitude, land_area_fraction, nc_vars = ga.nc_load(ncfile, ["altitude_of_0_degree_isotherm", "altitude_of_isoTprimW_equal_0"], time_period=hour_range)
In [5]:
print("From {0} to {1}".format(times[0], times[-1]))
len(times)
#plt.imshow(np.flipud(nc_vars['altitude_of_0_degree_isotherm'][6, :, :]))
Out[5]:
In [0]:
We use the parameters "altitude_of_0_degree_isotherm" and "altitude_of_isoTprimW_equal_0" from MEPS_extracted. Under dry conditions we use altitude_of_0_degree_isotherm and for timing we use the period with the highest values. With precipitation we use altitude_of_isoTprimW_equal_0 and the period with the highest amount of precipitation.
Logic: '''python If regional_precipitation < 2 (mm/døgn): use "altitude_of_0_degree_isotherm" else: find period of day with most intense precip (e.g. 06-12) use "altitude_of_isoTprimW_equal_0" for that period return freezing level '''
In [6]:
#fl_max = np.amax(nc_vars['altitude_of_0_degree_isotherm'][0:6,:,:], axis=0)
fl_max = np.amax(nc_vars['altitude_of_isoTprimW_equal_0'][0:24,:,:], axis=0)
fl_max
Out[6]:
In [7]:
plt.imshow(fl_max)
Out[7]:
In [44]:
# Load region mask - only for data on 1km xgeo-grid
region_id = 3031
region_mask, y_min, y_max, x_min, x_max = load_region(region_id)
print(y_max-y_min, x_max-x_min)
In [45]:
print(np.unique(region_mask))
plt.imshow(region_mask)
Out[45]:
In [46]:
var_name = 'altitude_of_isoTprimW_equal_0'
for t in range(len(times)):
t2 = t+hour_range[0]
t_str = times[t]
_fl = clip_region(np.flipud(nc_vars[var_name][t,:,:]), region_mask, t2, y_min, y_max, x_min, x_max)
plt.imshow(_fl, vmin=0, vmax=500, cmap='magma')
plt.axis('off')
plt.text(5, 5, "{0}: {1}".format(region_id, t_str), bbox=dict(facecolor='white', edgecolor='white', alpha=1.0))
cbar = plt.colorbar()
cbar.ax.set_ylabel(var_name)
_png = './img/fl_{1:02}.png'.format(region_id, t)
plt.savefig(_png, bbox_inches='tight')
plt.clf()
Use _makegif.py in folder img to generate a gif animation.
In [47]:
t_index = 0
fl_region = clip_region(np.flipud(fl_max), region_mask, t_index, y_min, y_max, x_min, x_max)
#print(np.unique(fl_3034))
In [48]:
#fl_3034.masked_where(fl_region<=0.0)
print(np.count_nonzero(np.isnan(fl_region)))
print(np.unique(fl_region))
plt.imshow(fl_region)
plt.colorbar()
Out[48]:
In [49]:
print("Mean\t: ", np.nanmean(fl_region.flatten()))
for p in [0,5,25,50,75, 80, 85,90, 95,100]:
print(p, "\t: ", np.nanpercentile(fl_region.flatten(), p))
In [50]:
fl_region_flat = fl_region[~np.isnan(fl_region)].data.flatten()
sns.distplot(fl_region_flat)
Out[50]:
In [51]:
nc_file2 = r"\\hdata\grid\metdata\prognosis\meps\det\archive\2019\meps_det_pp_1km_20190307T00Z.nc"
times, altitude, land_area_fraction, nc_vars2 = ga.nc_load(nc_file2, ["precipitation_amount"], time_period=hour_range)
In [58]:
precip_sum = np.sum(nc_vars2['precipitation_amount'][0:6,:,:], axis=0)
In [61]:
t_index = 0
precip_sum_region = clip_region(np.flipud(precip_sum), region_mask, t_index, y_min, y_max, x_min, x_max)
In [62]:
print(np.count_nonzero(np.isnan(precip_sum_region)))
print(np.unique(precip_sum_region))
plt.imshow(precip_sum_region)
plt.colorbar()
Out[62]:
In [65]:
# Mask where the precipitation during the day exceeds a given value.
psr_mask = np.where(precip_sum_region >= 5., 1, np.nan)
In [66]:
plt.imshow(psr_mask)
Out[66]:
In [67]:
fl_region_wet = fl_region * psr_mask
print("Mean\t: ", np.nanmean(fl_region_wet.flatten()))
for p in [0,5,25,50,75, 80, 85,90, 95,100]:
print(p, "\t: ", np.nanpercentile(fl_region_wet.flatten(), p))
In [ ]: