import numpy as np import matplotlib.pyplot as plt import matplotlib.colors as colors import cartopy.crs as ccrs from datetime import datetime from dateutil.relativedelta import relativedelta from utils import addCbar from utils_stations import loadStations, getStationsIniEnd from utils_cities import cities imgsDir = './imgs/paper/' #=============================================================================== # Load stations # #stations = loadStations('./stations/stationsData_20161201_to_20161202.pkl') #stations = loadStations('./stations/stationsData_20170101_to_20171231.pkl') #stations = loadStations('./stations/stationsData_20220101_to_20220831.pkl') #stations = loadStations('./stations/stationsData_20170101_to_20220831.pkl') stations = loadStations('./stations/stationsData_20170101_to_20221231.pkl') print('TODO: check whether I can get hourly measurements that are not averaged') #=============================================================================== # Data # ini_date, end_date = getStationsIniEnd(stations) years = list(np.arange(ini_date.year,(end_date+relativedelta(years=1)).year)) n_years = len(years) month_labels = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] month_labels_short = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'June', 'July', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] #=============================================================================== # Recompute some stuff for other time periods # ini_analysis = datetime(2017, 1,1) end_analysis = datetime(2023, 1,1) for nr, stn in stations.items(): idxs, times = stn.times_within(ini_analysis, end_analysis) # temperature bias stn_temps, t_idxs = stn.data('temp', times=times) if t_idxs != []: stn_sim_temps, _ = stn.data('sim_temp', times=list(np.array(times)[t_idxs])) # stn.mean_data['mean_diff'] = (stn_sim_temps - stn_temps).mean() stn.mean_data['mean_std'] = (stn_sim_temps - stn_temps).std() stn.mean_data['mean_corrected_diff'] = (stn_sim_temps + stn.t_corr - stn_temps).mean() stn.mean_data['mean_corrected_std'] = (stn_sim_temps + stn.t_corr - stn_temps).std() else: stn.mean_data['mean_diff'] = -99 stn.mean_data['mean_corrected_diff'] = -99 # precipitation bias stn_precs, t_idxs = stn.data('prec', times=times) if t_idxs != []: stn_sim_precs, _ = stn.data('sim_prec', times=list(np.array(times)[t_idxs])) # stn.mean_data['mean_prec_diff'] = (stn_sim_precs - stn_precs).mean() stn.mean_data['mean_prec_std'] = (stn_sim_precs - stn_precs).std() else: stn.mean_data['mean_prec_diff'] = -99 # rel humidity bias stn_rhums, t_idxs = stn.data('rhum', times=times) if t_idxs != []: stn_sim_rhums, _ = stn.data('sim_rhum', times=list(np.array(times)[t_idxs])) # stn.mean_data['mean_rhum_diff'] = (stn_sim_rhums - stn_rhums).mean() stn.mean_data['mean_rhum_std'] = (stn_sim_rhums - stn_rhums).std() else: stn.mean_data['mean_rhum_diff'] = -99 #=============================================================================== # Convert to arrays # nr_stn = len(stations) stn_STN = np.zeros(nr_stn) stn_lon = np.zeros(nr_stn) stn_lat = np.zeros(nr_stn) stn_bias = np.zeros(nr_stn) stn_hDelta = np.zeros(nr_stn) n_stn = 0 for nr, stn in stations.items(): if (stn.mean_data['mean_diff'] != -99): hDelta = stn.z_delta+2 # add the 2m for the temperature value if abs(hDelta) < 100: stn_STN [n_stn] = stn.STN stn_lon [n_stn] = stn.lon stn_lat [n_stn] = stn.lat stn_bias [n_stn] = stn.mean_data['mean_corrected_diff'] stn_hDelta[n_stn] = hDelta n_stn += 1 # trim stn_lon = stn_lon[0:n_stn] stn_lat = stn_lat[0:n_stn] stn_bias = stn_bias[0:n_stn] stn_hDelta = stn_hDelta[0:n_stn] stn_rhum_bias = np.zeros(nr_stn) n_stn = 0 for nr, stn in stations.items(): if (stn.mean_data['mean_rhum_diff'] != -99): hDelta = stn.z_delta+2 # add the 2m for the temperature value if abs(hDelta) < 100: stn_rhum_bias[n_stn] = stn.mean_data['mean_rhum_diff'] n_stn += 1 # trim stn_rhum_bias = stn_rhum_bias[0:n_stn] #=============================================================================== # Figure # average bias of simulations compared to weather stations # fig = plt.figure(1); plt.clf(); plt.show(block=False) ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) ax.set_title('Mean bias compared to weather stations {0:s} to {1:s}'.format(ini_analysis.strftime('%Y.%m.%d'),end_analysis.strftime('%Y.%m.%d'))) norm = colors.TwoSlopeNorm(vmin=stn_bias.min(), vcenter=0, vmax=stn_bias.max()) sc = ax.scatter(stn_lon, stn_lat, c=stn_bias, s=35, norm=norm, cmap='seismic', edgecolors='black') # plot city-rural stations differently idxs = [city['station_nr'] for city in cities] sc = ax.scatter([stn.lon for stn in [stations[id] for id in idxs]], \ [stn.lat for stn in [stations[id] for id in idxs]], \ c=[stn.mean_data['mean_diff'] for stn in [stations[id] for id in idxs]], \ marker='*', s=256, norm=norm, cmap='seismic', edgecolors='black') # add labels to all points for istn in range(n_stn): #label = '{0:d}: {1:d}'.format(int(stn_STN[istn]), int(stn_hDelta[istn])) label = '{0:d}'.format(int(stn_hDelta[istn])) plt.text(stn_lon[istn], stn_lat[istn], label, va='bottom', ha='right', fontsize='small') addCbar(fig, ax, sc, r'$<\Delta T>$ [K]') plt.savefig(imgsDir+'stations_map_{0:s}to{1:s}.png'.format(ini_analysis.strftime('%Y%m%d'),end_analysis.strftime('%Y%m%d'))) # #=============================================================================== # # Figure (included in the one above with labels) # # weather station height compared to grid cell height # # # fig = plt.figure(2); plt.clf(); plt.show(block=False) # ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) # ax.set_title('Weather station height - grid cell height') # #norm = colors.TwoSlopeNorm(vmin=stn_hDelta.min(), vcenter=0, vmax=stn_hDelta.max()) # #sc = ax.scatter(stn_lon, stn_lat, c=stn_hDelta, s=35, norm=norm, cmap='seismic') # norm = colors.TwoSlopeNorm(vmin=-50, vcenter=0, vmax=+50) # #sc = ax.scatter(stn_lon, stn_lat, c=stn_hDelta, s=35, norm=norm, cmap='seismic') # # divide markers # good_delta=10 # idxs = list(np.where(np.abs(stn_hDelta) <= good_delta)[0]) # #sc = ax.scatter(stn_lon, stn_lat, c=stn_hDelta, s=35, norm=norm, cmap='seismic') # sc = ax.scatter(stn_lon[idxs], stn_lat[idxs], c=stn_hDelta[idxs], marker='*', s=35, norm=norm, cmap='seismic') # addCbar(fig, ax, sc, r'$\Delta z$ [m]', extend='both') # #=============================================================================== # # Figure # # Time series at stations # # # data_label='rhum' # ini_analysis = datetime(2017, 1,1) # end_analysis = datetime(2022, 1,1) # fig = plt.figure(3); plt.clf(); plt.show(block=False) # fig, ax = plt.subplots(len(cities), 1, num=fig.number, sharex=True, sharey=False) # for icity, city in enumerate(cities): # stn = stations[city['station_nr']] # idxs, times = stn.times_within(ini_analysis, end_analysis) # # stn_data, t_idxs = stn.data(data_label, times=times) # t_array = np.array(times)[t_idxs] # stn_sim_data, _ = stn.data('sim_'+data_label, times=list(t_array)) # # if data_label == 'temp': # stn_data = stn_data + stn.t_corr # # ax[icity].plot(t_array, stn_data, label='Measurement') # ax[icity].plot(t_array, stn_sim_data, label='Simulated') # ax[icity].set_ylabel(city['name']) # ax[-1].legend(fontsize='small') # ax[-1].set_xlim([ini_analysis, end_analysis]) # plt.draw() # plt.savefig(imgsDir+'stations_timeHistory_{0:s}to{1:s}_{2:s}.png'.format(ini_analysis.strftime('%Y%m%d'),end_analysis.strftime('%Y%m%d'),data_label)) # #=============================================================================== # # Figure # # Time series at stations yearly COMPOSITE # # # ini_analysis = datetime(2017, 1,1) # end_analysis = datetime(2018, 1,1) # times = [ini_analysis + relativedelta(hours=x) for x in range((end_analysis-ini_analysis).days*24)] # nyears=5 # fig = plt.figure(31); plt.clf(); plt.show(block=False) # fig, ax = plt.subplots(len(cities), 1, num=fig.number, sharex=True, sharey=False) # for icity, city in enumerate(cities): # stn = stations[city['station_nr']] # # temps = np.zeros(len(times)) # sim_temps = np.zeros(len(times)) # counts = np.zeros(len(times)) # for y in range(nyears): # for itime, time in enumerate(times): # try: # temps [itime] += stn.time_data[time+relativedelta(years=y)]['temp'] # sim_temps[itime] += stn.time_data[time+relativedelta(years=y)]['sim_temp'] # counts [itime] += 1 # except: # print('not found') # # temps = temps / counts # sim_temps = sim_temps / counts # # ax[icity].plot(times, temps, label='Measurement') # ax[icity].plot(times, sim_temps+stn.t_corr, label='Simulated') # ax[icity].set_ylabel(city['name']) # ax[-1].legend(fontsize='small') # ax[-1].set_xlim([ini_analysis, end_analysis]) # plt.draw() # plt.savefig(imgsDir+'stations_timeHistory_composite.png') #=============================================================================== # Figure # Time series at stations SEASON DIEL COMPOSITE # for icity, city in enumerate(cities): stn = stations[city['station_nr']] #print('{0:2d} {1:10s}: {2:.2f}C | deltaZ: {3:.2f}'.format(icity, city['name'], stn.mean_data['mean_corrected_diff'], stn.z_delta)) print('{0:2d} {1:10s}: {2:5.2f}C | {3:10.2e} | {4:7.2f} | deltaZ: {5:7.2f}'.format( \ icity, city['name'], stn.mean_data['mean_corrected_diff'], stn.mean_data['mean_prec_diff'], stn.mean_data['mean_rhum_diff'], stn.z_delta \ )) #s_cities = [cities[7], cities[6]] # St. Gallen, Lucerne (fucked up approx lat/lon) s_cities = [cities[0], cities[8]] # Zurich, Lugano ini_analysis = datetime(2017, 1,1) end_analysis = datetime(2022, 1,1) fig = plt.figure(32, figsize=(8, 4.55)); plt.clf(); plt.show(block=False) fig, ax = plt.subplots(len(s_cities), 4, num=fig.number, sharex=True, sharey='row') for icity, city in enumerate(s_cities): stn = stations[city['station_nr']] times = np.arange(24) temps = np.zeros((4, 24)) sim_temps = np.zeros((4, 24)) counts = np.zeros((4, 24)) v_temps = np.zeros((4, 24)) v_sim_temps = np.zeros((4, 24)) idxs, _ = stn.times_within(ini_analysis, end_analysis) a_times = stn.times() a_temps = stn.temps() a_sim_temps = stn.sim_temps() # Mean for idx in idxs: month = a_times[idx].month if (month >=12) | (month<=2): season = 0 elif (month > 2) & (month<=5): season = 1 elif (month > 5) & (month<=8): season = 2 else: season = 3 hour = a_times[idx].hour temps [season, hour] += a_temps [idx] sim_temps[season, hour] += a_sim_temps[idx] + stn.t_corr counts [season, hour] += 1 temps = temps / counts sim_temps = sim_temps / counts # Variance for idx in idxs: month = a_times[idx].month if (month >=12) | (month<=2): season = 0 elif (month > 2) & (month<=5): season = 1 elif (month > 5) & (month<=8): season = 2 else: season = 3 hour = a_times[idx].hour v_temps [season, hour] += (a_temps [idx] - temps[season,hour])**2 v_sim_temps[season, hour] += (a_sim_temps[idx] + stn.t_corr - sim_temps[season,hour])**2 v_temps = ( v_temps / (counts-1))**0.5 v_sim_temps = (v_sim_temps / (counts-1))**0.5 for season in range(4): ax[icity,season].plot(times, temps[season,:], label='Measurement') ax[icity,season].plot(times, sim_temps[season,:], label='Simulated') ax[icity,season].fill_between(times, temps[season,:]- v_temps[season,:], temps[season,:]+ v_temps[season,:],alpha=0.5) ax[icity,season].fill_between(times, sim_temps[season,:]-v_sim_temps[season,:], sim_temps[season,:]+v_sim_temps[season,:],alpha=0.5) ax[icity,season].set_xlim([0,23]) ax[icity,season].set_xticks([0, 6, 12, 18]) ax[icity,0].set_ylabel(city['name']) # save city[ 'times'] = times city[ 'temps'] = temps city[ 'sim_temps'] = sim_temps city[ 'v_temps'] = v_temps city['v_sim_temps'] = v_sim_temps ax[-1,0].set_xlabel('Winter [DJF]') ax[-1,1].set_xlabel('Spring [MAM]') ax[-1,2].set_xlabel('Summer [JJA]') ax[-1,3].set_xlabel('Autumn [SON]') ax[-1,-1].legend(fontsize='small') plt.draw() plt.savefig(imgsDir+'stations_timeHistory_seasonComposite.png') plt.savefig(imgsDir+'stations_timeHistory_seasonComposite.pdf') # Try to quantify phase #import scipy.fftpack # Number of samplepoints N = 24 # sample spacing T = 1.0 / 24 fig = plt.figure(332); plt.clf(); plt.show(block=False) fig, ax = plt.subplots(2,4, num=fig.number, sharex=True,sharey=True) for icity, city in enumerate(s_cities): for season in range(4): x = city['times'] y0 = city[ 'temps'][season,:] y1 = city['sim_temps'][season,:] delta_mea = y0.max() - y0.min() delta_sim = y1.max() - y1.min() diff = delta_mea-delta_sim print('delta mea = {0:.2f} | delta sim = {1:.2f} | delta = {2:.2f} | rerr = {3:.2f}'.format(delta_mea, delta_sim, diff, diff/delta_mea)) xf = np.linspace(0.0, 1.0/(2.0*T), N//2) #y0f = scipy.fftpack.fft(y0) #y1f = scipy.fftpack.fft(y1) #ax[icity,season].plot(xf, 2.0/N * np.abs(y0f[:N//2])) #ax[icity,season].plot(xf, 2.0/N * np.abs(y1f[:N//2])) plt.draw() #=============================================================================== # Figure # Wind rose # from windrose import WindroseAxes from matplotlib import cm stn = stations[71] nrows, ncols = 1, 2 nbins=8 nsector=32 normed=True lw=2 fig = plt.figure(333, figsize=(13,6)); plt.clf(); plt.show(block=False) ax = fig.add_subplot(nrows, ncols, 1, projection="windrose") ax.contour( stn.data( 'winD'), stn.data( 'winM'), normed=normed, bins=np.arange(0, nbins, 1), nsector=nsector, cmap=cm.autumn, lw=lw) ax.set_legend() ax = fig.add_subplot(nrows, ncols, 2, projection="windrose") ax.contour((180+stn.data('sim_winD'))%360, stn.data('sim_winM'), normed=normed, bins=np.arange(0, nbins, 1), nsector=nsector, cmap=cm.winter, lw=lw) ax.set_legend() plt.draw() #=============================================================================== # Figure # Time series of daily max UHI # fig = plt.figure(4, figsize=(12,8)); plt.clf(); plt.show(block=False) fig, ax = plt.subplots(3, 1, num=fig.number, sharex=True, sharey=False) for icity, city in enumerate(cities): stn = stations[city['station_nr']] all_times = np.array(stn.times()) all_sims = np.array(stn.sim_temps()) all_cities = np.array(stn.city_temps()) days = [] uhis = [] p0 = p1 = 0 while p1 < len(all_times)-24: p0 = p1 try: p1 = np.argwhere(all_times >= all_times[p0]+relativedelta(days=1))[0][0] except: p1 = len(all_times)+1 days.append( all_times[p0] ) uhis.append( np.max( all_cities[p0:p1] - all_sims[p0:p1] ) ) days = np.array(days) uhis = np.array(uhis) ax[0].plot(days, uhis, label=city['name']) # monthly mean value mo_days = [] mo_uhis = [] p0 = p1 = 0 while p1 < len(days): p0 = p1 try: p1 = np.argwhere(days >= days[p0]+relativedelta(months=1))[0][0] except: p1 = len(days)+1 mo_days.append( days[p0] ) mo_uhis.append( np.mean( uhis[p0:p1] ) ) ax[1].plot(mo_days, mo_uhis, label=city['name']) # yearly mean value ye_days = [] ye_uhis = [] p0 = p1 = 0 while p1 < len(days): p0 = p1 try: p1 = np.argwhere(days >= days[p0]+relativedelta(years=1))[0][0] except: p1 = len(days)+1 ye_days.append( days[p0] ) ye_uhis.append( np.mean( uhis[p0:p1] ) ) ax[2].plot(ye_days, ye_uhis, label=city['name']) ax[1].grid('on') ax[2].grid('on') ax[-1].legend(fontsize='small') plt.draw() plt.savefig(imgsDir+'uhi_timeSeries_atStations') #=============================================================================== # Figure # diel cicle of temperature # hours = np.arange(24) for year in years: fig = plt.figure(11, figsize=(8,13)); plt.clf(); plt.show(block=False) fig, ax = plt.subplots(12, len(cities), num=fig.number, sharex=True, sharey='row') fig.suptitle('Diel cycle of 2m temperature for the year {0:d}'.format(year)) for month in range(12): ax[month, 0].set_ylabel(month_labels_short[month]) s_date = datetime(year,month+1,1) e_date = s_date+relativedelta(months=1) for icity, city in enumerate(cities): diel_mea = np.zeros(24) diel_sim = np.zeros(24) diel_cit = np.zeros(24) diel_simC = np.zeros(24) counters = np.zeros(24) stn = stations[city['station_nr']] all_times = np.array(stn.times()) all_temps = np.array(stn.temps()) all_sims = np.array(stn.sim_temps()) all_cities = np.array(stn.city_temps()) all_simsC = np.array(stn.sim_temps()+stn.t_corr) idxs = list(np.argwhere( (all_times>=s_date) & (all_times 0: diel_mea [hour] = diel_mea [hour] / counters[hour] diel_sim [hour] = diel_sim [hour] / counters[hour] diel_cit [hour] = diel_cit [hour] / counters[hour] diel_simC[hour] = diel_simC[hour] / counters[hour] ax[month, icity].plot(hours, diel_mea, label='Rural measurement') ax[month, icity].plot(hours, diel_sim, label='Rural simulated') ax[month, icity].plot(hours, diel_cit, label='City simulated') ax[month, icity].plot(hours, diel_simC, '--', label='Rural simulated C') if month == 0: ax[month, icity].set_title(city['name']) elif month == 11: ax[month, icity].set_xticks([0, 6, 12, 18]) ax[-1,-1].legend(fontsize='small') plt.draw() plt.savefig(imgsDir+'diel_cycle_{0:d}.png'.format(year)) #=============================================================================== # Figure # diel cycle of UHI # hours = np.arange(24) for year in years: fig = plt.figure(12, figsize=(8,13)); plt.clf(); plt.show(block=False) fig, ax = plt.subplots(12, len(cities), num=fig.number, sharex=True, sharey='row') fig.suptitle('Diel cycle of UHI for the year {0:d}'.format(year)) for month in range(12): ax[month, 0].set_ylabel(month_labels_short[month]) s_date = datetime(year,month+1,1) e_date = s_date+relativedelta(months=1) for icity, city in enumerate(cities): diel_uhi = np.zeros(24) counters = np.zeros(24) stn = stations[city['station_nr']] all_times = np.array(stn.times()) all_sims = np.array(stn.sim_temps()) all_cities = np.array(stn.city_temps()) idxs = list(np.argwhere( (all_times>=s_date) & (all_times 0: diel_uhi[hour] = diel_uhi[hour] / counters[hour] ax[month, icity].plot(hours, diel_uhi, label='UHI') if month == 0: ax[month, icity].set_title(city['name']) elif month == 11: ax[month, icity].set_xticks([0, 6, 12, 18]) #ax[-1,-1].legend(fontsize='small') plt.draw() plt.savefig(imgsDir+'diel_uhi_{0:d}.png'.format(year)) #=============================================================================== # Figure # Tropical nights # fig = plt.figure(13); plt.clf(); plt.show(block=False) fig, ax = plt.subplots(2, 1, num=fig.number, sharex=True, sharey=True) #fig.suptitle('Number of tropical nights per year') scale=0.095 for icity, city in enumerate(cities): ax[0].bar(np.array(years) - scale*(len(cities)-1)/2 + scale*icity, city['sta_rural_nights'], width = scale, label=city['name']) ax[1].bar(np.array(years) - scale*(len(cities)-1)/2 + scale*icity, city['sta_city_nights'], width = scale, label=city['name']) ax[0].set_ylabel('Rural areas') ax[1].set_ylabel('Cities') ax[1].set_xticks(years) ax[0].grid(which='major', axis='y') ax[1].grid(which='major', axis='y') ax[0].set_axisbelow(True) ax[1].set_axisbelow(True) # Shrink the axes for iax in range(len(ax)): box = ax[iax].get_position() ax[iax].set_position([box.x0, box.y0, box.width * 0.9, box.height]) # Put a legend to the right of the current axis ax[0].legend(loc='center left', bbox_to_anchor=(1, 0), fontsize='small') #ax[0].legend(fontsize='small') plt.draw() plt.savefig(imgsDir+'tropical_nights_atStations')