Homec4science
Diffusion solar_potential (aw_solar) Browse

solarPOT/libs/visualise_spacetime.py
eee6149cc313aw_solar

Last Change
Raw File

visualise_spacetime.py
View Options

import numpy as np
import pandas as pd
import xarray as xr
import os
import time
import util
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.animation as manimation
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.widgets import Slider, Button, RadioButtons
import calendar
# COLORMAP FOR SOLAR RADIATION VISUALISATION
colorsRGB = [(47,111,191),
(51,148,190),
(65,195,153),
(105,219,77),
(183,233,81),
(249,239,88),
(250,201,53),
(246,135,26),
(220,82,7),
(193,31,4)]
colors = []
for row in colorsRGB:
colors.append(tuple([i/255.0 for i in row]))
cm = LinearSegmentedColormap.from_list('radiation', colors, N=500)
def check_df(arr, cols):
if arr.shape[1] in (cols, cols-1):
return arr
else: return arr.T
def plot_map_year( filename = None, data = None,
h0 = None, m0 = None,
interactive = True, animation = True, savefile = None,
var = 'SIS', vmin = 0, vmax = 1000,
ctext = None, colormap = None,
scale = 1):
# scale font sizes and figure shapes
BASE_FONT = 15*scale
TIMELABEL_FONT = 20*scale
FIG_WIDTH = 15*scale
FIG_HEIGHT = 7.5*scale
# check/load data
if filename is None and data is None:
print('Please provide some input data')
print('EXIT')
return
if data is None:
data = xr.open_dataset(filename)
# create data frames
data_frames = []
months = []
hours = []
for month in range(12):
for hour in range(len(data[var].hour.values)):
if np.sum(data[var].isel(date = month, hour = hour)>0.0001):
data_frames.append(data[var].isel(date = month, hour = hour))
months.append(int(month+1))
if month < 4 or month > 10:
hours.append(int(hour+4))
else:
hours.append(int(hour+5))
# xarray dataset with frames
df = pd.DataFrame(data = np.vstack([months, hours]).T, columns = ['month', 'hour']).reset_index().set_index(['month', 'hour'])
frames = df.to_xarray().fillna(0)
frames.index.values = frames.index.values.astype('int')
frames = frames.index
if interactive:
if m0 is None:
m0 = 6
if h0 is None:
h0 = 13
elif animation:
if m0 is None:
m0 = months[0]
if h0 is None:
h0 = hours[0]
else:
if m0 is None or h0 is None:
print('Please pass month (m0) and/or hour(h0) to be plotted\nEXIT')
return
f0 = int(frames.sel(month = m0, hour = h0))
# define update functions
def update_slider(val):
mon = smonth.val
hr = shour.val
f = int(frames.sel(month = mon, hour = hr))
update(f)
def update(frame_number):
Z = data_frames[frame_number].values[:-1, :-1]
Zm = np.ma.masked_invalid(Z)
p.set_array(check_df(Zm,nx).ravel())
txt.set_text("%s\n\n%02d:00h" %(calendar.month_name[months[frame_number]], hours[frame_number]))
# PLOTTING
if ctext is None:
ctext = 'Solar radiation (in $W/m^2$)'
if colormap is None:
colormap = cm
matplotlib.rcParams.update({'font.size': BASE_FONT})
fig, ax = plt.subplots(figsize = (FIG_WIDTH,FIG_HEIGHT))
ax.axis('equal')
ax.axis('off')
nx = len(data.x.values)
ny = len(data.y.values)
X, Y = np.meshgrid(data.x.values, data.y.values,indexing='xy')
Z = data_frames[f0].values[:-1, :-1]
Zm = np.ma.masked_invalid(Z)
# plt.title("Day: %d, hour: %d" %(day+1, hour+3), horizontalalignment='right')
txt = ax.text(X[0,0], Y[-1,-1]-30000, ("%s\n\n%02d:00h" %(calendar.month_name[m0], h0)), color='black', style='italic', fontsize = TIMELABEL_FONT)
p = ax.pcolormesh(X,Y,check_df(Zm,nx), vmin = vmin, vmax = vmax)
p.set_cmap(colormap)
cbar = plt.colorbar(p)
cbar.set_label(ctext, labelpad=0)
labels = [item.get_text() for item in cbar.ax.get_yticklabels()]
labels[-1] = '> ' + labels[-1]
cbar.ax.set_yticklabels(labels)
plt.show()
if interactive:
fig.subplots_adjust(bottom=0.18, top=0.93)
axcolor = 'whitesmoke'
axmonth = plt.axes([0.2, 0.05, 0.5, 0.03], facecolor=axcolor)
axhour = plt.axes([0.2, 0.1, 0.5, 0.03], facecolor=axcolor)
smonth = Slider(axmonth, 'Month', 1, 12, valinit=6, valfmt = '%d', valstep = 1, color='silver')
smonth.vline.set_visible(False)
shour = Slider(axhour, 'Hour', 5, 21, valinit=13, valfmt = '%d', valstep = 1, color='silver')
shour.vline.set_visible(False)
smonth.on_changed(update_slider)
shour.on_changed(update_slider)
plt.show()
elif animation:
print("Generating video with %d frames" %len(data_frames))
my_animation = manimation.FuncAnimation(fig, update, frames = len(data_frames), interval=200, repeat = False)
if savefile is None:
plt.draw()
else:
timer = util.Timer()
my_animation.save(savefile)
timer.stop()
else:
if savefile is not None: plt.savefig(savefile)
def plot_map_month(filename = None, data = None,
month = 6, h0 = None,
interactive = True, animation = True, savefile = None,
var = 'SIS',
vmin = 0, vmax = 1000,
ctext = None, colormap = None,
scale = 1,
title = None):
# scale font sizes and figure shapes
BASE_FONT = 15*scale
TIMELABEL_FONT = 20*scale
FIG_WIDTH = 15*scale
FIG_HEIGHT = 7.5*scale
# check/load data
if filename is None and data is None:
print('Please provide some input data')
print('EXIT')
return
if data is None:
data = xr.open_dataset(filename)
# create data frames
data_frames = []
hours = []
for hour in range(len(data[var].hour.values)):
if np.sum(data[var].isel(hour = hour)>0.0001):
data_frames.append(data[var].isel(hour = hour))
if month < 4 or month > 10:
hours.append(int(hour+4))
else:
hours.append(int(hour+5))
df = pd.DataFrame(data = hours, columns = ['hour']).reset_index().set_index(['hour'])
frames = df.to_xarray()
frames = frames.index
if interactive:
if h0 is None:
h0 = 13
elif animation:
if h0 is None:
h0 = index[0,1]
else:
if h0 is None:
print('Please pass hour(h0) to be plotted\nEXIT')
return
f0 = int(frames.sel(hour = h0))
# define update functions
def update_slider(val):
hr = shour.val
f = int(frames.sel(hour = hr))
update(f)
def update(frame_number):
Z = data_frames[frame_number].values[:-1, :-1]
Zm = np.ma.masked_invalid(Z)
p.set_array(check_df(Zm,nx).ravel())
txt.set_text("%s\n\n%02d:00h" %(calendar.month_name[month], hours[frame_number]))
# PLOTTING
if ctext is None:
ctext = 'Solar radiation (in $W/m^2$)'
if colormap is None:
colormap = cm
matplotlib.rcParams.update({'font.size': BASE_FONT})
fig, ax = plt.subplots(figsize = (FIG_WIDTH,FIG_HEIGHT))
ax.axis('equal')
ax.axis('off')
nx = len(data.x.values)
ny = len(data.y.values)
X, Y = np.meshgrid(data.x.values, data.y.values)
Z = data_frames[f0].values[:-1, :-1]
Zm = np.ma.masked_invalid(Z)
# plt.title("Day: %d, hour: %d" %(day+1, hour+3), horizontalalignment='right')
txt = ax.text(X[0,0], Y[-1,-1]-30000, ("%s\n\n%02d:00h" %(calendar.month_name[month], h0)), color='black', style='italic', fontsize = TIMELABEL_FONT)
p = ax.pcolormesh(X,Y,check_df(Zm,nx), vmin = vmin, vmax = vmax)
p.set_cmap(colormap)
cbar = plt.colorbar(p)
cbar.set_label(ctext, labelpad=0)
labels = [item.get_text() for item in cbar.ax.get_yticklabels()]
labels[-1] = '> ' + labels[-1]
cbar.ax.set_yticklabels(labels)
if title is not None:
fig.subplots_adjust(top=0.83)
plt.title(title, position = (0.5,1.1), fontsize = int(TIMELABEL_FONT * 1.2))
plt.show()
if interactive:
fig.subplots_adjust(bottom=0.18, top=0.93)
axcolor = 'whitesmoke'
axhour = plt.axes([0.2, 0.1, 0.5, 0.03], facecolor=axcolor)
shour = Slider(axhour, 'Hour', 5, 21, valinit=13, valfmt = '%d', valstep = 1, color='silver')
shour.vline.set_visible(False)
shour.on_changed(update_slider)
plt.show()
elif animation:
print("Generating video with %d frames" %len(data_frames))
my_animation = manimation.FuncAnimation(fig, update, frames = len(data_frames), interval=200, repeat = False)
if savefile is None:
plt.draw()
else:
timer = util.Timer()
my_animation.save(savefile)
timer.stop()
else:
if savefile is not None:
plt.savefig(savefile, dpi=300)
print('Saved figure as %s with resolution 300 dpi' %savefile)
def plot_unc_month(path, month = 6, h0 = None, interactive = True, animation = True, savefile = None, var = 'SIS', vmin = 0, vmax = 1000):
pred = xr.open_dataset(os.path.join(path, 'prediction.nc'))
model_unc = xr.open_dataset(os.path.join(path, 'model_unc.nc'))
data_unc = xr.open_dataset(os.path.join(path, 'data_unc.nc'))
total_unc = xr.open_dataset(os.path.join(path, 'total_unc.nc'))
# create data frames
data_frames_pred = []
data_frames_model = []
data_frames_data = []
data_frames_unc = []
hours = []
for hour in range(17):
if np.sum(data[var].isel(hour = hour)>0.0001):
data_frames_pred.append(pred[var].isel(hour = hour))
data_frames_model.append(model_unc[var].isel(hour = hour))
data_frames_data.append(data_unc[var].isel(hour = hour))
data_frames_unc.append(total_unc[var].isel(hour = hour))
if month < 4 or month > 10:
hours.append(int(hour+4))
else:
hours.append(int(hour+5))
df = pd.DataFrame(data = hours, columns = ['hour']).reset_index().set_index(['hour'])
frames = df.to_xarray()
frames = frames.index
if interactive:
if h0 is None:
h0 = 13
elif animation:
if h0 is None:
h0 = index[0,1]
else:
if h0 is None:
print('Please pass hour(h0) to be plotted\nEXIT')
return
f0 = int(frames.sel(hour = h0))
# define update functions
def update_slider(val):
hr = shour.val
f = int(frames.sel(hour = hr))
update(f)
def update(frame_number):
Z_p = data_frames_pred[frame_number].values[:-1, :-1]
Z_m = data_frames_model[frame_number].values[:-1, :-1]
Z_d = data_frames_data[frame_number].values[:-1, :-1]
Z_u = data_frames_unc[frame_number].values[:-1, :-1]
Zm_p = np.ma.masked_invalid(Z_p)
Zm_m = np.ma.masked_invalid(Z_m)
Zm_d = np.ma.masked_invalid(Z_d)
Zm_u = np.ma.masked_invalid(Z_u)
p.set_array(check_df(Zm,nx).ravel())
txt.set_text("%s\n\n%02d:00h" %(calendar.month_name[month], hours[frame_number]))
# PLOTTING
if ctext is None:
ctext = 'Solar radiation (in $W/m^2$)'
if colormap is None:
colormap = cm
matplotlib.rcParams.update({'font.size': 22})
fig, ax = plt.subplots(figsize = (20,10))
ax.axis('equal')
ax.axis('off')
nx = len(data.x.values)
ny = len(data.y.values)
X, Y = np.meshgrid(data.x.values, data.y.values)
Z = data_frames[f0].values[:-1, :-1]
Zm = np.ma.masked_invalid(Z)
# plt.title("Day: %d, hour: %d" %(day+1, hour+3), horizontalalignment='right')
txt = ax.text(X[0,0], Y[-1,-1]-30000, ("%s\n\n%02d:00h" %(calendar.month_name[month], h0)), color='black', style='italic', fontsize = 28)
p = ax.pcolormesh(X,Y,check_df(Zm,nx), vmin = vmin, vmax = vmax)
p.set_cmap(colormap)
cbar = plt.colorbar(p)
cbar.set_label(ctext, labelpad=0)
labels = [item.get_text() for item in cbar.ax.get_yticklabels()]
labels[-1] = '> ' + labels[-1]
cbar.ax.set_yticklabels(labels)
plt.show()
if interactive:
fig.subplots_adjust(bottom=0.18, top=0.93)
axcolor = 'whitesmoke'
axhour = plt.axes([0.2, 0.1, 0.5, 0.03], facecolor=axcolor)
shour = Slider(axhour, 'Hour', 5, 21, valinit=13, valfmt = '%d', valstep = 1, color='silver')
shour.vline.set_visible(False)
shour.on_changed(update_slider)
plt.show()
elif animation:
print("Generating video with %d frames" %len(data_frames))
my_animation = manimation.FuncAnimation(fig, update, frames = len(data_frames), interval=200, repeat = False)
if savefile is None:
plt.draw()
else:
timer = util.Timer()
my_animation.save('RF_mmh_200ms.mp4')
timer.stop()
else:
if savefile is not None: plt.savefig(savefile)
c4science · Help