# Licensed to the Apache Software Foundation (ASF) under one | |
# or more contributor license agreements. See the NOTICE file | |
# distributed with this work for additional information | |
# regarding copyright ownership. The ASF licenses this file | |
# to you under the Apache License, Version 2.0 (the | |
# "License"); you may not use this file except in compliance | |
# with the License. You may obtain a copy of the License at | |
# | |
# http://www.apache.org/licenses/LICENSE-2.0 | |
# | |
# Unless required by applicable law or agreed to in writing, | |
# software distributed under the License is distributed on an | |
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY | |
# KIND, either express or implied. See the License for the | |
# specific language governing permissions and limitations | |
# under the License. | |
#Apache OCW lib immports | |
import ocw.dataset as ds | |
import ocw.data_source.local as local | |
import ocw.plotter as plotter | |
import ocw.utils as utils | |
from ocw.evaluation import Evaluation | |
import ocw.metrics as metrics | |
# Python libraries | |
import numpy as np | |
import numpy.ma as ma | |
import matplotlib.pyplot as plt | |
from mpl_toolkits.basemap import Basemap | |
from matplotlib import rcParams | |
from matplotlib.patches import Polygon | |
import string | |
def Map_plot_bias_of_multiyear_climatology(obs_dataset, obs_name, model_datasets, model_names, | |
file_name, row, column, map_projection=None): | |
'''Draw maps of observed multi-year climatology and biases of models"''' | |
# calculate climatology of observation data | |
obs_clim = utils.calc_temporal_mean(obs_dataset) | |
# determine the metrics | |
map_of_bias = metrics.TemporalMeanBias() | |
# create the Evaluation object | |
bias_evaluation = Evaluation(obs_dataset, # Reference dataset for the evaluation | |
model_datasets, # list of target datasets for the evaluation | |
[map_of_bias, map_of_bias]) | |
# run the evaluation (bias calculation) | |
bias_evaluation.run() | |
rcm_bias = bias_evaluation.results[0] | |
fig = plt.figure() | |
lat_min = obs_dataset.lats.min() | |
lat_max = obs_dataset.lats.max() | |
lon_min = obs_dataset.lons.min() | |
lon_max = obs_dataset.lons.max() | |
string_list = list(string.ascii_lowercase) | |
nmodels = len(model_datasets) | |
row, column = plotter._best_grid_shape(nmodels + 1, (row, column)) | |
ax = fig.add_subplot(row,column,1) | |
if map_projection == 'npstere': | |
m = Basemap(ax=ax, projection ='npstere', boundinglat=lat_min, lon_0=0, | |
resolution = 'l', fix_aspect=True) | |
else: | |
m = Basemap(ax=ax, projection ='cyl', llcrnrlat = lat_min, urcrnrlat = lat_max, | |
llcrnrlon = lon_min, urcrnrlon = lon_max, resolution = 'l', fix_aspect=True) | |
if obs_dataset.lons.ndim == 1 and obs_dataset.lats.ndim == 1: | |
lons, lats = np.meshgrid(obs_dataset.lons, obs_dataset.lats) | |
if obs_dataset.lons.ndim == 2 and obs_dataset.lats.ndim == 2: | |
lons = obs_dataset.lons | |
lats = obs_dataset.lats | |
x,y = m(lons, lats) | |
m.drawcoastlines(linewidth=1) | |
m.drawcountries(linewidth=1) | |
m.drawstates(linewidth=0.5, color='w') | |
max = m.contourf(x,y,obs_clim,levels = plotter._nice_intervals(obs_dataset.values, 10), extend='both',cmap='rainbow') | |
ax.annotate('(a) \n' + obs_name,xy=(lon_min, lat_min)) | |
cax = fig.add_axes([0.02, 1.-float(1./row)+1./row*0.25, 0.01, 1./row*0.5]) | |
plt.colorbar(max, cax = cax) | |
clevs = plotter._nice_intervals(rcm_bias, 11) | |
for imodel in np.arange(nmodels): | |
ax = fig.add_subplot(row, column,2+imodel) | |
if map_projection == 'npstere': | |
m = Basemap(ax=ax, projection ='npstere', boundinglat=lat_min, lon_0=0, | |
resolution = 'l', fix_aspect=True) | |
else: | |
m = Basemap(ax=ax, projection ='cyl', llcrnrlat = lat_min, urcrnrlat = lat_max, | |
llcrnrlon = lon_min, urcrnrlon = lon_max, resolution = 'l', fix_aspect=True) | |
m.drawcoastlines(linewidth=1) | |
m.drawcountries(linewidth=1) | |
m.drawstates(linewidth=0.5, color='w') | |
max = m.contourf(x,y,rcm_bias[imodel,:],levels = clevs, extend='both', cmap='RdBu_r') | |
ax.annotate('('+string_list[imodel+1]+') \n '+model_names[imodel],xy=(lon_min, lat_min)) | |
cax = fig.add_axes([0.91, 0.5, 0.015, 0.4]) | |
plt.colorbar(max, cax = cax) | |
plt.subplots_adjust(hspace=0.01,wspace=0.05) | |
fig.savefig(file_name,dpi=600,bbox_inches='tight') | |
def Taylor_diagram_spatial_pattern_of_multiyear_climatology(obs_dataset, obs_name, model_datasets, model_names, | |
file_name): | |
# calculate climatological mean fields | |
obs_clim_dataset = ds.Dataset(obs_dataset.lats, obs_dataset.lons, obs_dataset.times, utils.calc_temporal_mean(obs_dataset)) | |
model_clim_datasets = [] | |
for dataset in model_datasets: | |
model_clim_datasets.append(ds.Dataset(dataset.lats, dataset.lons, dataset.times, utils.calc_temporal_mean(dataset))) | |
# Metrics (spatial standard deviation and pattern correlation) | |
# determine the metrics | |
taylor_diagram = metrics.SpatialPatternTaylorDiagram() | |
# create the Evaluation object | |
taylor_evaluation = Evaluation(obs_clim_dataset, # Climatological mean of reference dataset for the evaluation | |
model_clim_datasets, # list of climatological means from model datasets for the evaluation | |
[taylor_diagram]) | |
# run the evaluation (bias calculation) | |
taylor_evaluation.run() | |
taylor_data = taylor_evaluation.results[0] | |
plotter.draw_taylor_diagram(taylor_data, model_names, obs_name, file_name, pos='upper right',frameon=False) | |
def Time_series_subregion(obs_subregion_mean, obs_name, model_subregion_mean, model_names, seasonal_cycle, | |
file_name, row, column, x_tick=['']): | |
nmodel, nt, nregion = model_subregion_mean.shape | |
if seasonal_cycle: | |
obs_data = ma.mean(obs_subregion_mean.reshape([1,nt/12,12,nregion]), axis=1) | |
model_data = ma.mean(model_subregion_mean.reshape([nmodel,nt/12,12,nregion]), axis=1) | |
nt = 12 | |
else: | |
obs_data = obs_subregion_mean | |
model_data = model_subregion_mean | |
x_axis = np.arange(nt) | |
x_tick_values = x_axis | |
fig = plt.figure() | |
rcParams['xtick.labelsize'] = 6 | |
rcParams['ytick.labelsize'] = 6 | |
for iregion in np.arange(nregion): | |
ax = fig.add_subplot(row, column, iregion+1) | |
x_tick_labels = [''] | |
if iregion+1 > column*(row-1): | |
x_tick_labels = x_tick | |
else: | |
x_tick_labels=[''] | |
ax.plot(x_axis, obs_data[0, :, iregion], color='r', lw=2, label=obs_name) | |
for imodel in np.arange(nmodel): | |
ax.plot(x_axis, model_data[imodel, :, iregion], lw=0.5, label = model_names[imodel]) | |
ax.set_xlim([-0.5,nt-0.5]) | |
ax.set_xticks(x_tick_values) | |
ax.set_xticklabels(x_tick_labels) | |
ax.set_title('Region %02d' % (iregion+1), fontsize=8) | |
ax.legend(bbox_to_anchor=(-0.2, row/2), loc='center' , prop={'size':7}, frameon=False) | |
fig.subplots_adjust(hspace=0.7, wspace=0.5) | |
fig.savefig(file_name, dpi=600, bbox_inches='tight') | |
def Portrait_diagram_subregion(obs_subregion_mean, obs_name, model_subregion_mean, model_names, seasonal_cycle, | |
file_name, normalize=True): | |
nmodel, nt, nregion = model_subregion_mean.shape | |
if seasonal_cycle: | |
obs_data = ma.mean(obs_subregion_mean.reshape([1,nt/12,12,nregion]), axis=1) | |
model_data = ma.mean(model_subregion_mean.reshape([nmodel,nt/12,12,nregion]), axis=1) | |
nt = 12 | |
else: | |
obs_data = obs_subregion_mean | |
model_data = model_subregion_mean | |
subregion_metrics = ma.zeros([4, nregion, nmodel]) | |
for imodel in np.arange(nmodel): | |
for iregion in np.arange(nregion): | |
# First metric: bias | |
subregion_metrics[0, iregion, imodel] = metrics.calc_bias(model_data[imodel, :, iregion], obs_data[0, :, iregion], average_over_time = True) | |
# Second metric: standard deviation | |
subregion_metrics[1, iregion, imodel] = metrics.calc_stddev_ratio(model_data[imodel, :, iregion], obs_data[0, :, iregion]) | |
# Third metric: RMSE | |
subregion_metrics[2, iregion, imodel] = metrics.calc_rmse(model_data[imodel, :, iregion], obs_data[0, :, iregion]) | |
# Fourth metric: correlation | |
subregion_metrics[3, iregion, imodel] = metrics.calc_correlation(model_data[imodel, :, iregion], obs_data[0, :, iregion]) | |
if normalize: | |
for iregion in np.arange(nregion): | |
subregion_metrics[0, iregion, : ] = subregion_metrics[0, iregion, : ]/ma.std(obs_data[0, :, iregion])*100. | |
subregion_metrics[1, iregion, : ] = subregion_metrics[1, iregion, : ]*100. | |
subregion_metrics[2, iregion, : ] = subregion_metrics[2, iregion, : ]/ma.std(obs_data[0, :, iregion])*100. | |
region_names = ['R%02d' % i for i in np.arange(nregion)+1] | |
for imetric, metric in enumerate(['bias','std','RMSE','corr']): | |
plotter.draw_portrait_diagram(subregion_metrics[imetric, :, :], region_names, model_names, file_name+'_'+metric, | |
xlabel='model',ylabel='region') | |
def Map_plot_subregion(subregions, ref_dataset, directory): | |
lons, lats = np.meshgrid(ref_dataset.lons, ref_dataset.lats) | |
fig = plt.figure() | |
ax = fig.add_subplot(111) | |
m = Basemap(ax=ax, projection='cyl',llcrnrlat = lats.min(), urcrnrlat = lats.max(), | |
llcrnrlon = lons.min(), urcrnrlon = lons.max(), resolution = 'l') | |
m.drawcoastlines(linewidth=0.75) | |
m.drawcountries(linewidth=0.75) | |
m.etopo() | |
x, y = m(lons, lats) | |
#subregion_array = ma.masked_equal(subregion_array, 0) | |
#max=m.contourf(x, y, subregion_array, alpha=0.7, cmap='Accent') | |
for subregion in subregions: | |
draw_screen_poly(subregion[1], m, 'w') | |
plt.annotate(subregion[0],xy=(0.5*(subregion[1][2]+subregion[1][3]), 0.5*(subregion[1][0]+subregion[1][1])), ha='center',va='center', fontsize=8) | |
fig.savefig(directory+'map_subregion', bbox_inches='tight') | |
def draw_screen_poly(boundary_array, m, linecolor='k'): | |
''' Draw a polygon on a map | |
:param boundary_array: [lat_north, lat_south, lon_east, lon_west] | |
:param m : Basemap object | |
''' | |
lats = [boundary_array[0], boundary_array[0], boundary_array[1], boundary_array[1]] | |
lons = [boundary_array[3], boundary_array[2], boundary_array[2], boundary_array[3]] | |
x, y = m( lons, lats ) | |
xy = zip(x,y) | |
poly = Polygon( xy, facecolor='none',edgecolor=linecolor ) | |
plt.gca().add_patch(poly) |