diff --git a/FIGSEXAMPLES/demo_slides_destine-obsall-apps.pdf b/FIGSEXAMPLES/demo_slides_destine-obsall-apps.pdf new file mode 100644 index 0000000000000000000000000000000000000000..47ae73eec8e78b54c9b713d834ed7f696ec88ed3 Binary files /dev/null and b/FIGSEXAMPLES/demo_slides_destine-obsall-apps.pdf differ diff --git a/FIGSEXAMPLES/demo_video_destine-obsall-apps.mp4 b/FIGSEXAMPLES/demo_video_destine-obsall-apps.mp4 new file mode 100644 index 0000000000000000000000000000000000000000..24128ba2a9343ce28327f38173c082a7293d00d2 Binary files /dev/null and b/FIGSEXAMPLES/demo_video_destine-obsall-apps.mp4 differ diff --git a/FIGSEXAMPLES/destine_obs_radsound_figs.jpg b/FIGSEXAMPLES/destine_obs_radsound_figs.jpg index b6330ad151e98685c03b424610c796c13e2e0652..6956adddbc1d4941a485d707ec4ae75de8fea170 100644 Binary files a/FIGSEXAMPLES/destine_obs_radsound_figs.jpg and b/FIGSEXAMPLES/destine_obs_radsound_figs.jpg differ diff --git a/FIGSEXAMPLES/destine_obs_synop_figs.jpg b/FIGSEXAMPLES/destine_obs_synop_figs.jpg index 7db4166b71c12c4ea9a3f2aaa58098707ae3372c..a453472944614c45c143b48ba305d51cf2117a7e 100644 Binary files a/FIGSEXAMPLES/destine_obs_synop_figs.jpg and b/FIGSEXAMPLES/destine_obs_synop_figs.jpg differ diff --git a/FIGSEXAMPLES/front_demo_destine-obsall-apps.jpg b/FIGSEXAMPLES/front_demo_destine-obsall-apps.jpg new file mode 100644 index 0000000000000000000000000000000000000000..355362d68ea3c1a41dd8f35c4ef95de674124b65 Binary files /dev/null and b/FIGSEXAMPLES/front_demo_destine-obsall-apps.jpg differ diff --git a/RADSOUND/STATRS/fortran-programs/rank_histogram_summary_statistics.f95 b/RADSOUND/STATRS/fortran-programs/rank_histogram_summary_statistics.f95 index 46a8a608af26ddb7e8e1cd32fc327b9aac536c28..6b59109b7123815681cbda21bd0d03029f6e780d 100644 --- a/RADSOUND/STATRS/fortran-programs/rank_histogram_summary_statistics.f95 +++ b/RADSOUND/STATRS/fortran-programs/rank_histogram_summary_statistics.f95 @@ -11,7 +11,7 @@ ! that the input file (INFILE) only includes data for these ! hours. ! -! Jouni Räisänen, University of Helsinki, August 2023 +! Jouni Räisänen, University of Helsinki, October 2023 ! !------------------------------------------------------------------------- ! @@ -250,6 +250,11 @@ max_freq_p=0 max_freq_q=0 do i=1,nbins_p + do hour=hour1,hour2,hour_step + if(frequency_p(i,hour+1).gt.max_freq_p/nbins_p)then + max_freq_p=frequency_p(i,hour+1)*nbins_p + endif + enddo if(frequency_p(i,nhour+1).gt.max_freq_p/nbins_p)then max_freq_p=frequency_p(i,nhour+1)*nbins_p endif @@ -257,6 +262,11 @@ max_freq_p=1.02*max_freq_p do i=1,nquant+1 + do hour=hour1,hour2,hour_step + if(frequency_q(i,hour+1).gt.max_freq_q/(nquant+1))then + max_freq_q=frequency_q(i,hour+1)*(nquant+1) + endif + enddo if(frequency_q(i,nhour+1).gt.max_freq_q/(nquant+1))then max_freq_q=frequency_q(i,nhour+1)*(nquant+1) endif @@ -272,10 +282,8 @@ open(unit=1,form='formatted',file=text_output) write(1,'(A20,I3)')'Number of stations: ',n_station write(1,'(A24,I3)')'Number of p-value bins: ',nbins_p - if(grads_out)then - write(1,'(A23,F6.3)')'Axis_limit_p_for_GrADS ',max_freq_p - write(1,'(A23,F6.3)')'Axis_limit_q_for_GrADS ',max_freq_q - endif + write(1,'(A23,F6.3)')'Axis_limit_p ',max_freq_p + write(1,'(A23,F6.3)')'Axis_limit_q ',max_freq_q do hour=hour1,hour2,hour_step write(1,'(A22,I2.2,A6,F5.3)')'Frequency, p <= 0.001 ',hour,' UTC: ',& frequency_p01(hour+1) @@ -319,12 +327,16 @@ ! Write the file giving the p-values for individual stations ! open(unit=2,form='formatted',status='unknown',file=outfile_p_values) - headerline=' station@hdr:integer longitude@hdr:real latitude@hdr:real p_value@body:real' + headerline=' station@hdr:integer longitude@hdr:real latitude@hdr:real hour@hdr:integer p_value@body:real' write(2,*)trim(headerline) do i=1,n_station - write(2,'(7X,I6,3F16.6)')& ! this far: just fwrite the station ordinal number, in case - ! the station code was a character string - station(i),longitude(i),latitude(i),p_value(nhour+1,i) + do hour=hour1,hour2,hour_step + write(2,'(7X,I6,2F16.6,I3,F16.6)')& ! this far: just fwrite the station ordinal number, in case + ! the station code was a character string + station(i),longitude(i),latitude(i),hour, p_value(hour+1,i) + enddo + write(2,'(7X,I6,2F16.6,I3,F16.6)')& + station(i),longitude(i),latitude(i),24,p_value(nhour+1,i) enddo close(2) ! diff --git a/RADSOUND/STATRS/python/plot_p_values_map.py b/RADSOUND/STATRS/python/plot_p_values_map.py deleted file mode 100644 index 12442d9428dadd2cdfc1a0f540c5afbdd58a7b48..0000000000000000000000000000000000000000 --- a/RADSOUND/STATRS/python/plot_p_values_map.py +++ /dev/null @@ -1,88 +0,0 @@ -#!/usr/bin/env python - -''' -Written by Madeleine Ekblom, 2023-04-06 - -Based on Jouni Räisänen's script for plotting p values on the map of Finland - -Plots p values on a map restricted to an area limited by lons (19, 32) and lats (59, 71) - -Example: -$ python3 plot_p_values_map.py p_values_as_text_file -$ python3 plot_p_values_map.py p-values - -p_values_as_text_file is a text file containing stationID, longitude, latitude, and p value with one header row -''' - -import numpy as np -import matplotlib.pyplot as plt -import sys -import cartopy.crs as ccrs -from matplotlib.lines import Line2D - -def read_data(p_values_file): - ''' - Description: Reads in file containing p values at different stations - Input: file containing p values with stationID (sid), longitude (lon), latitude (lat), and p value (p). Header = 1 row. - Output: structure numpy array with sid, lon, lat, p - ''' - - p_values_data = np.loadtxt(p_values_file, skiprows=1, dtype={'names': ('sid', 'lon', 'lat', 'p'), 'formats':('i4', 'f4', 'f4', 'f4')}) - - return p_values_data - - -def plot_p_values(p_values): - ''' - Description: plot p values on a map restricted to an area limited by lons (19, 32) and lats (59, 71) - Input: numpy array - Output: saves a figure in the running directory - ''' - - lon_min=-180 - lon_max=180 - lat_min=-90 - lat_max=90 - - colors = ['darkviolet', 'blue', 'skyblue', 'lime', 'yellow', 'orange', 'red', 'grey'] - limits = [0,0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,1] - - fig = plt.figure(figsize=(12,4.5)) - ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree()) - ax.coastlines() - ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree()) - ax.gridlines(draw_labels=True) - - # basic plot, not changing colors of p: - #ax.scatter(p_values['lon'], p_values['lat'], c=p_values['p']) - - #loop over limits: - legend_elements = [] - for n in range(8): - p_ind = ((p_values['p'] > limits[n]) & (p_values['p'] <= limits[n+1])) - ax.scatter(p_values['lon'][p_ind], p_values['lat'][p_ind], c=colors[n]) - legend_elements.append(Line2D([0], [0], marker='o', color='w',markerfacecolor=colors[n], label=str(limits[n])+'-'+str(limits[n+1]))) - - # legend table with customized values - - #print(legend_elements) - - ax.legend(handles=legend_elements, loc='center left', bbox_to_anchor=(1.1, 0.5)) - ax.set_title('P-values for T2-T850 quantiles in ' + timestring, fontsize=20) - - plt.savefig('p_values.png', dpi=300) -# plt.show() - - -def main(p_values_file): - - p_values_data = read_data(p_values_file) - - plot_p_values(p_values_data) - - -if __name__=='__main__': - p_values_file = sys.argv[1] - timestring = sys.argv[2] - - main(p_values_file) diff --git a/RADSOUND/STATRS/python/plot_p_values_map_rank_hist_temp.py b/RADSOUND/STATRS/python/plot_p_values_map_rank_hist_temp.py new file mode 100644 index 0000000000000000000000000000000000000000..11c8b55a20754458b0888bd2e527a963d60acf04 --- /dev/null +++ b/RADSOUND/STATRS/python/plot_p_values_map_rank_hist_temp.py @@ -0,0 +1,139 @@ +#!/usr/bin/env python + +''' +Written by Jouni Räisänen, 2023-10-07 +Updates by Madeleine Ekblom, 2023-11-29 +- move legend into lower left corner of plot +- figure title added, subplot titles changed + +Plots p values from quantile bin rank histogram analysis to global maps (scale 0-1, only low values indicate significance). +P-values are plotted separately for 00 UTC and 12 UTC, to take into account varying data availability. + +Example: +$ python3 plot_p_values_map_rank_hist_temp.py p_values_as_text_file_00 p_values_as_text_file_12 timestring figname + +p_values_as_text_file_XX are text files containing stationID, longitude, latitude, and p value with one header row +timestring indicates the simulation period used. It is only used in the map headers. +''' + +import numpy as np +import matplotlib.pyplot as plt +import sys +import cartopy.crs as ccrs +from matplotlib.lines import Line2D + +def read_data(p_values_file): + ''' + Description: Reads in file containing p values at different stations + Input: file containing p values with stationID (sid), longitude (lon), latitude (lat), and p value (p). Header = 1 row. + Output: structure numpy array with sid, lon, lat, p + ''' + + p_values_data = np.loadtxt(p_values_file, skiprows=1, dtype={'names': ('sid', 'lon', 'lat', 'p'), 'formats':('S3', 'f4', 'f4', 'f4')}) + + return p_values_data + + +def plot_p_values_00_12(p_values_00, p_values_12): + ''' + Description: plot p values at 00 and 12 UTC on global maps + Input: two numpy arrays + Output: saves a figure in the running directory + ''' + + # define lat and lon for global map + lon_min=-180 + lon_max=180 + lat_min=-90 + lat_max=90 + + # define colour and limits for the dots + colors = ['darkviolet', 'blue', 'skyblue', 'lime', 'yellow', 'orange', 'red', 'grey'] + limits = [0,0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,1] + + # set up figure + fig = plt.figure(figsize=(12,11)) + + # plotting details: + fs1=20 # fontsize for title + fs2=12 # fontsize for labels, subplot titles + fs3=10 # fontsize for figure explanation + + ## 00 UTC ## + # add subplot with map, for 00 UTC data + ax = fig.add_subplot(2,1,1,projection=ccrs.PlateCarree()) + ax.coastlines() + ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree()) + ax.gridlines(draw_labels=True) + + #loop over limits: + legend_elements = [] + for n in range(8): + # find indices for p-values + p_mod = 0.5+1.0001*(p_values_00['p']-0.5) + p_ind = ((p_mod > limits[n]) & (p_mod <= limits[n+1])) + # scatter plot with dots coloured based on p-values + ax.scatter(p_values_00['lon'][p_ind], p_values_00['lat'][p_ind], c=colors[n]) + # add to legend_elements: colour and limit + legend_elements.append(Line2D([0], [0], marker='o', color='w',markerfacecolor=colors[n], label=str(limits[n])+'-'+str(limits[n+1]))) + + # add legend to lower left corner of map + ax.legend(handles=legend_elements, loc='lower left') + + # add subplot title + ax.set_title('00 UTC', fontsize=fs2, fontweight='bold') + + ## 12 UTC ## + # add subplot with map, for 12 UTC data + ax = fig.add_subplot(2,1,2,projection=ccrs.PlateCarree()) + ax.coastlines() + ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree()) + ax.gridlines(draw_labels=True) + + #loop over limits: + legend_elements = [] + for n in range(8): + # find indices for p-values + p_mod = 0.5+1.0001*(p_values_12['p']-0.5) + p_ind = ((p_mod > limits[n]) & (p_mod <= limits[n+1])) + # scatter plot with dots coloured based on p-values + ax.scatter(p_values_12['lon'][p_ind], p_values_12['lat'][p_ind], c=colors[n]) + # add to legend elements: colour and limit + legend_elements.append(Line2D([0], [0], marker='o', color='w',markerfacecolor=colors[n], label=str(limits[n])+'-'+str(limits[n+1]))) + + # add legend to lower left corner of map + ax.legend(handles=legend_elements, loc='lower left') + + # add subplot title + ax.set_title('12 UTC', fontsize=fs2, fontweight='bold') + + # explanation for figure: + text_expl ='p-value: the probability that the MSD statistics for the quantile rank \nhistogram would be exceeded by pure chance when comparing \na perfect model simulation with observations.' + + # add text box with explanation to the bottom of the figure + ax.text(0.2, -0.25, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3, transform=ax.transAxes) + + # add figure title: + fig.suptitle('RADSOUND stations: p-values for $\Delta$T = 2t - t850\nquantiles in ' + timestring, fontsize=fs1) #, fontweight='bold') + + # save figure + plt.savefig(figname, dpi=300) + + +def main(p_values_file_00, p_values_file_12): + ''' + Main: (1) reads in p values at 00 UTC and 12 UTC (2) plots p-values on a global map + ''' + + p_values_data_00 = read_data(p_values_file_00) + p_values_data_12 = read_data(p_values_file_12) + + plot_p_values_00_12(p_values_data_00, p_values_data_12) + +if __name__=='__main__': + p_values_file_00 = sys.argv[1] + p_values_file_12 = sys.argv[2] + timestring = sys.argv[3] + figname = sys.argv[4] + + main(p_values_file_00, p_values_file_12) diff --git a/RADSOUND/STATRS/python/plot_quantiles_rankhist.py b/RADSOUND/STATRS/python/plot_quantiles_rankhist.py deleted file mode 100644 index 5478aa0e4bb1a98fd8c82059a62a9d98b8b8b063..0000000000000000000000000000000000000000 --- a/RADSOUND/STATRS/python/plot_quantiles_rankhist.py +++ /dev/null @@ -1,193 +0,0 @@ -#!/usr/bin/env python - -''' -Written by Madeleine Ekblom, 2023-05-10 - -Based on Jouni Räisänen's script for plotting quantiles/time series data (00, 06, 12, 18 UTC) and rank histogram (24UTC=combination of 00,06,12,and 18 UTC data) - -Example: -$ python3 plot_quantiles_rankhist.py example.cfg - -example.cfg is a text file containing: -[data] -variable=39 -station_id=102019 -longitude=23.576000 -latitude=68.602997 -year_beg=2010 -year_end=2012 -month_beg=01 -month_end=12 - -[input] -sim_data=time_series_102019_2010-2012.txt -quantile_sel=quantiles_102019.txt -rank_hists=rank_histogram_102019_2010-2012.txt - -[output] -fig_name=standard_plot_102019_2010-2012_python.png - -''' - -import numpy as np -import matplotlib.pyplot as plt -import sys -import configparser -import datetime -from matplotlib.dates import drange, MonthLocator, DateFormatter - -config = configparser.ConfigParser() - -def plot_data(quantiles_data, rank_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, savefig=False): - ''' - quantiles data give the lines and the time series data give the dots - ''' - # set up date array consisting of all days in beg year: - days = np.arange(1,366) -# day1 = datetime.date(year_beg, 1, 1) -# day2 = datetime.date(year_beg+1, 1, 1) - day1 = datetime.date(2001, 1, 1) - day2 = datetime.date(2002, 1, 1) - delta = datetime.timedelta(days=1) - dates = drange(day1, day2, delta) - - # arrays for quantile names, axes titles, hours to plot - q_names = ['q01', 'q10', 'q25', 'q50', 'q75', 'q90', 'q99'] - sub_titles = ['00 UTC', '06 UTC', '12 UTC', '18 UTC'] - hours = [0, 6, 12, 18] - - # calculate number of years: - nyears = year_end - year_beg + 1 - - # set up figure - fig = plt.figure(figsize=(10,10), layout='constrained') - spec = fig.add_gridspec(3,2) - fig.suptitle(fig_title, fontsize=20) - - # counter - c = 0 - - # plot quantiles/time series for times 00, 06, 12, 18: - for i in range(2): - for j in range(2): - # set up axis: title, xaxis - ax = fig.add_subplot(spec[i,j]) - ax.set_title(sub_titles[c]) - ax.set_xlim(dates[0], dates[-1]) - ax.xaxis.set_major_locator(MonthLocator()) - ax.xaxis.set_major_formatter(DateFormatter('%b')) - - # quantile data: - # find quantile data hour = 0, 6, 12, 18 - qh_ind = (quantiles_data['hour'][:] == hours[c]) - for q in q_names: - q_data = quantiles_data[q][qh_ind] - ax.plot(dates, q_data[:365], 'k-') # --> change q only contains one year of data!! - # plot time series data: - # find time series where hour = 0, 6, 12, 18 - th_ind = (time_series_data['hour'][:] == hours[c]) - t_data = time_series_data['value'][th_ind] # all years - for n in range(nyears): - ax.scatter(dates, t_data[n*365:(n+1)*365], marker='.') - c = c + 1 - - # plot rank histogram data: - ax2 = fig.add_subplot(spec[2,:]) - ax2.set_title('Rank histogram\n' + 'MSD:' + str(rank_data[4,4]) + ' p-value:' + str(rank_data[4,5])) - ax2.bar(np.arange(100), rank_data[4,6:]*100) - ax2.set_xlim(0,100) - ax2.axhline(y=1) - - if savefig: - plt.savefig(fig_name, dpi=300) - -# plt.show( ) - - -def read_quantiles_data(quantile_sel_file): - ''' - Reads quantile data from text file and returns a structured numpy array - Input: text file - Output:a structured numpy array: sid, lon, lat, day_of_year, hour, q01, q10, q25, q50, q75, q90, q99 - ''' - # Header line contains: - #station@hdr:integer longitude@hdr:real latitude@hdr:real day_of_year@hdr:integer hour@hdr:integer q01@body:real q10@body:real q25@body:real q50@body:real q75@body:real q90@body:real q99@body:real - - quantile_data = np.loadtxt(quantile_sel_file, skiprows=1, dtype={'names':('sid', 'lon', 'lat', 'day_of_year', 'hour', 'q01', 'q10', 'q25', 'q50', 'q75', 'q90', 'q99'), 'formats': ('i4', 'f4', 'f4', 'i4', 'i4', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4')}) - - return quantile_data - -def read_rank_data(rank_hists_file): - ''' - Reads rank histogram data binned into a 100 bins: 0-1%, 1-2%, ..., 99-100% - Input: text file - Output: numpy array (unstructured) first columns contain station id, longitude, latitude, hour, msd, p_value, and then follows the bins f00 = 0-1%, ..., f99=99-100% - ''' - - #Header line contains: - # station@hdr:integer longitude@hdr:real latitude@hdr:real hour@hdr:integer msd@body:real p_value@body:real f00@body:real f01@body:real f02@body:real f03@body:real f04@body:real f05@body:real f06@body:real f07@body:real f08@body:real f09@body:real f10@body:real f11@body:real f12@body:real f13@body:real f14@body:real f15@body:real f16@body:real f17@body:real f18@body:real f19@body:real f20@body:real f21@body:real f22@body:real f23@body:real f24@body:real f25@body:real f26@body:real f27@body:real f28@body:real f29@body:real f30@body:real f31@body:real f32@body:real f33@body:real f34@body:real f35@body:real f36@body:real f37@body:real f38@body:real f39@body:real f40@body:real f41@body:real f42@body:real f43@body:real f44@body:real f45@body:real f46@body:real f47@body:real f48@body:real f49@body:real f50@body:real f51@body:real f52@body:real f53@body:real f54@body:real f55@body:real f56@body:real f57@body:real f58@body:real f59@body:real f60@body:real f61@body:real f62@body:real f63@body:real f64@body:real f65@body:real f66@body:real f67@body:real f68@body:real f69@body:real f70@body:real f71@body:real f72@body:real f73@body:real f74@body:real f75@body:real f76@body:real f77@body:real f78@body:real f79@body:real f80@body:real f81@body:real f82@body:real f83@body:real f84@body:real f85@body:real f86@body:real f87@body:real f88@body:real f89@body:real f90@body:real f91@body:real f92@body:real f93@body:real f94@body:real f95@body:real f96@body:real f97@body:real f98@body:real f99@body:real - - rank_hists = np.loadtxt(rank_hists_file, skiprows=1) - - return rank_hists - -def read_time_series_data(sim_data_file): - ''' - Reads time series data - Input: text file - Output: structured numpy array: sid, lon, lat, year, day_of_year, hour, value - ''' - - # Header line contains - # station@hdr:integer longitude@hdr:real latitude@hdr:real year@hdr:integer day_of_year@hdr:integer hour@hdr:integer value@body:real - - time_series_data = np.loadtxt(sim_data_file, skiprows=1, dtype={'names': ('sid', 'lon', 'lat', 'year', 'day_of_year', 'hour', 'value'), 'formats': ('i4', 'f4', 'f4', 'i4', 'i4', 'i4', 'f4')}) - - return time_series_data - -def main(config_file): - ''' - Main: - reads config files - reads input data: time_series_data, quantiles_data, rank_hist_data - calls plotting function and saves figure with figure name given in config file - ''' - ## Read from config file ## - config.read(config_file) - - # data to be plotted - variable = config['data']['variable'] - station_id = config['data']['station_id'] - longitude = config['data']['longitude'] - latitude = config['data']['latitude'] - year_beg = int(config['data']['year_beg']) - year_end = int(config['data']['year_end']) - month_beg = config['data']['month_beg'] - month_end = config['data']['month_end'] - - # input files - sim_file = config['input']['sim_data'] - quantile_file = config['input']['quantile_sel'] - rank_hist_file = config['input']['rank_hists'] - - # output files - fig_name = config['output']['fig_name'] - - ## Read input data ## - time_series_data = read_time_series_data(sim_file) - quantiles_data = read_quantiles_data(quantile_file) - rank_hist_data = read_rank_data(rank_hist_file) - - ## Plot data ## - fig_title = 'Station: ' + station_id + ', ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end + '\nLat=' + latitude + ' Lon=' + longitude - plot_data(quantiles_data, rank_hist_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, True) - - -if __name__=='__main__': - if (len(sys.argv) < 2): - sys.exit("Error: config file must be added as 2nd argument") - elif (len(sys.argv) > 2): - sys.exit("Error: only add config file as argument") - config_file = sys.argv[1] - - main(config_file) diff --git a/RADSOUND/STATRS/python/plot_quantiles_rankhist_00_12.py b/RADSOUND/STATRS/python/plot_quantiles_rankhist_00_12.py index b9d632ebe3c70d7ad8a395acda612670efb97897..003391a2f03439444a95ff7246520c8c18b9f115 100644 --- a/RADSOUND/STATRS/python/plot_quantiles_rankhist_00_12.py +++ b/RADSOUND/STATRS/python/plot_quantiles_rankhist_00_12.py @@ -2,8 +2,12 @@ ''' Written by Madeleine Ekblom, 2023-05-10 +Updates 2023-12-01: +- visual changes of the plot: labels, legends, title, and colours updated +- text boxes added with rank histogram MSD and p-values +- added: check if data available, if not add text box "No data available" -Based on Jouni Räisänen's script for plotting quantiles/time series data (00, 06, 12, 18 UTC) and rank histogram (24UTC=combination of 00,06,12,and 18 UTC data) +Based on Jouni Räisänen's script for plotting quantiles/time series data for SYNOP data, modified to plot only 00, 12 UTC data and rank histogram at 00, 12 UTC Example: $ python3 plot_quantiles_rankhist.py example.cfg @@ -35,82 +39,212 @@ import sys import configparser import datetime from matplotlib.dates import drange, MonthLocator, DateFormatter +import matplotlib.lines as mlines +import matplotlib.patches as mpatches +import matplotlib.ticker as ticker config = configparser.ConfigParser() -def plot_data(quantiles_data, rank_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, savefig=False): +def plot_data(quantiles_data, rank_data, time_series_data, station_id, lon, lat, year_beg, year_end, month_beg, month_end, fig_name, savefig=False): ''' - quantiles data give the lines and the time series data give the dots + Plot data at 00UTC and 12UTC: quantiles data give the lines and the time series data give the dots. ''' # set up date array consisting of all days in beg year: days = np.arange(1,366) -# day1 = datetime.date(year_beg, 1, 1) -# day2 = datetime.date(year_beg+1, 1, 1) + # day1 = datetime.date(year_beg, 1, 1) + # day2 = datetime.date(year_beg+1, 1, 1) day1 = datetime.date(2001, 1, 1) day2 = datetime.date(2002, 1, 1) delta = datetime.timedelta(days=1) dates = drange(day1, day2, delta) + + + + # experiment id: now use 'test data' - to be replaced with real exp_id (as input argument) + exp_id='test data' # arrays for quantile names, axes titles, hours to plot q_names = ['q01', 'q10', 'q25', 'q50', 'q75', 'q90', 'q99'] - sub_titles = ['00 UTC', '06 UTC', '12 UTC', '18 UTC'] - hours = [0, 6, 12, 18] + sub_titles = ['00 UTC', '12 UTC'] + hours = [0, 12] # calculate number of years: nyears = year_end - year_beg + 1 - # set up figure + # set up figure with 4x2 subplots -> last two rows of subplots for legend and textboxes fig = plt.figure(figsize=(10,10), layout='constrained') - spec = fig.add_gridspec(3,2) + spec = fig.add_gridspec(4,2) + + # plotting details: + alpha= [0.5, 0.3, 0.1] # transparency for shaded area plotting + color='blue' # color for shading + linecolor='red' # color for solid line + scatter_size=5 # size for scattered dots + fs1=20 # fontsize for title + fs2=14 # fontsize for labels, subplot titles + fs3=12 # fontsize for legend + fs4=10 # fontize for figure textbox + ndigits=3 # round off to ndigits (e.g., lon, lat) + + # legend elemnents for quantile data, one legend for model and one legend for observations + legend_elements1 = [mlines.Line2D([],[], marker='.', linestyle='', color='k', label='Model data (' + exp_id + ')')] + legend_elements2= [mpatches.Patch(color=color, alpha=alpha[0], label='25-50%'), + mpatches.Patch(color=color, alpha=alpha[1], label='10-25%'), + mpatches.Patch(color=color, alpha=alpha[2], label='1-10%'), + mlines.Line2D([],[], color=linecolor, label='50%'), + mpatches.Patch(color=color, alpha=alpha[0], label='50-75%'), + mpatches.Patch(color=color, alpha=alpha[1], label='75-90%'), + mpatches.Patch(color=color, alpha=alpha[2], label='90-99%')] + + # text box for empty data: + text_no_data = "No data available" + + # round lon, lat to three digits: + lon = np.round(float(lon), ndigits) + lat = np.round(float(lat), ndigits) + + # define lon/lat text strings for title of the form lat: XX N/S, lat: XX E/W + if (lat > 0): + # if lat positive -> N + lat_text = ' (lat: ' + str(lat) + '$^\circ$N, ' + else: + # if negative -> S + lat_text = ' (lat ' + str(-lat) + '$^\circ$S, ' + + if (0 <= lon) and (lon <= 180): + # lon btw 0 an 180 -> E + lon_text = ' lon: ' + str(lon) + '$^\circ$E)' + else: + if lon < 0: + # lon < 0 --> - lon W + lon_text = ' lon: ' + str(-lon) + '$^\circ$W)' + else: + # 180 < lon < 360 --> - (lon - 360) W + lon = -(lon - 360) + lon_text = ' lon: ' + str(lon) + '$^\circ$W)' + + # define figure title: + fig_title = 'RADSOUND station: ' + station_id + ' ' + lat_text + lon_text + '\nModel data from ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end fig.suptitle(fig_title, fontsize=20) # counter c = 0 - # plot quantiles/time series for times 00, 06, 12, 18: - for i in range(1): - for j in range(2): - # set up axis: title, xaxis - ax = fig.add_subplot(spec[i,j]) - ax.set_title(sub_titles[c]) - ax.set_xlim(dates[0], dates[-1]) - ax.xaxis.set_major_locator(MonthLocator()) - ax.xaxis.set_major_formatter(DateFormatter('%b')) - - # quantile data: - # find quantile data hour = 0, 6, 12, 18 - qh_ind = (quantiles_data['hour'][:] == hours[c]) - for q in q_names: - q_data = quantiles_data[q][qh_ind] - ax.plot(dates, q_data[:365], 'k-') # --> change q only contains one year of data!! - # plot time series data: - # find time series where hour = 0, 6, 12, 18 - th_ind = (time_series_data['hour'][:] == hours[c]) - t_data = time_series_data['value'][th_ind] # all years + # plot quantiles/time series for times 00, 12. Loop from 0-1 to plot first row of subplots: + for i in range(2): + # set up axis: title, xaxis, x- and y-labels + ax = fig.add_subplot(spec[0,i]) + ax.set_title(sub_titles[c], fontsize=fs2) + ax.set_xlim(dates[0], dates[-1]) + ax.xaxis.set_major_locator(MonthLocator()) + ax.xaxis.set_major_formatter(DateFormatter('%b')) + ax.set_xlabel('Time', fontsize=fs2) + ax.set_ylabel('$\Delta$T = 2t - t850 ($^\circ$C)', fontsize=fs2) + + # quantile data: + # find quantile data hour = 0, 12 + qh_ind = (quantiles_data['hour'][:] == hours[c]) + + # get quantile data, qXX_data = quantile XX data at hour qh_ind + q01_data = quantiles_data[q_names[0]][qh_ind] + q10_data = quantiles_data[q_names[1]][qh_ind] + q25_data = quantiles_data[q_names[2]][qh_ind] + q50_data = quantiles_data[q_names[3]][qh_ind] + q75_data = quantiles_data[q_names[4]][qh_ind] + q90_data = quantiles_data[q_names[5]][qh_ind] + q99_data = quantiles_data[q_names[6]][qh_ind] + + # solid line at q50: + ax.plot(dates, q50_data, c=linecolor, ls='-') + + # shading will be the most opaque close to 50 and more transparent close to 1/99 + + # fill area btw q25-q50, q50-q75: alpha = alpha[0] + ax.fill_between(dates, q25_data[:], q50_data[:], color=color, alpha=alpha[0]) + ax.fill_between(dates, q50_data[:], q75_data[:], color=color, alpha=alpha[0]) + + # fill area btw q10-q25, q75-q90 + ax.fill_between(dates, q10_data[:], q25_data[:], color=color, alpha=alpha[1]) + ax.fill_between(dates, q75_data[:], q90_data[:], color=color, alpha=alpha[1]) + + # fill area btw q01-q10, q90-q99 + ax.fill_between(dates, q01_data[:], q10_data[:], color=color, alpha=alpha[2]) + ax.fill_between(dates, q90_data[:], q99_data[:], color=color, alpha=alpha[2]) + + # plot time series data: + # find time series where hour = 0, 12 + th_ind = (time_series_data['hour'][:] == hours[c]) + t_data = time_series_data['value'][th_ind] # all years + + # check if data available for given station and time: check that the length of array t_data > 0 + if (len(t_data) > 0): for n in range(nyears): - ax.scatter(dates, t_data[n*365:(n+1)*365], marker='.') - c = c + 2 - - # plot rank histogram data: - ax2 = fig.add_subplot(spec[1,:]) - ax2.set_title('Rank histogram\n' + 'MSD:' + str(rank_data[0,4]) + ' p-value:' + str(rank_data[0,5])) - ax2.bar(np.arange(100), rank_data[0,6:]*100) - ax2.set_xlim(0,100) - ax2.axhline(y=1) - - ax2 = fig.add_subplot(spec[2,:]) - ax2.set_title('Rank histogram\n' + 'MSD:' + str(rank_data[2,4]) + ' p-value:' + str(rank_data[2,5])) - ax2.bar(np.arange(100), rank_data[2,6:]*100) - ax2.set_xlim(0,100) - ax2.axhline(y=1) - + # plot data with black, small dots: + ax.scatter(dates, t_data[n*365:(n+1)*365], s=scatter_size, c='k', marker='.') + else: + ax.text(0.5, 0.5, text_box, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + + # increase counter by 1 + c = c + 1 + + # plot rank histogram data, loop from 0-1 to plot last row of subplots: + # rank_data[0,6:] contains 00UTC and rank_data[2,6:] contains 12UTC data --> plot rank_data[i*2, 6:] + for i in range(2): + ax2 = fig.add_subplot(spec[1,i]) + ax2.set_title('Rank histogram for station: ' + station_id + ' at ' + sub_titles[i], fontsize=fs2) + + # check if data available for given station and time: if MSD value = -999 000 then no data. MSD value = rank_data[i*2=0, 12 UTC data, 4] > -99999 + if (rank_data[i*2,4] > -99999): + ax2.bar(np.arange(100)+0.5, rank_data[i*2,6:]*100) # add + 0.5 to center bars + else: + ax2.text(0.5, 0.5, text_no_data, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + ax2.set_xlim(0,100) + ax2.axhline(y=1, c='k', ls='--') # horizontal black dashed line at y=0 + ax2.set_xlabel('Percentile', fontsize=fs2) + ax2.set_ylabel('Relative frequency', fontsize=fs2) + ax2.xaxis.set_major_locator(ticker.FixedLocator([10,25,50,75,90])) + + # third row of axis: + # axis for legend, position: 2x0 + ax3 = fig.add_subplot(spec[2,0]) + ax3.axis('off') + + # add two legends: model data, and observations (with two columns) + leg1 = ax3.legend(loc='upper center', handles=legend_elements1, fontsize=fs4) + leg2 = ax3.legend(loc='lower center', handles=legend_elements2, ncols=2, fontsize=fs4, title='Quantiles for observations', title_fontsize=fs4) + ax3.add_artist(leg1) + ax3.add_artist(leg2) + + # axis for rank histogram text, positition 2x1 + ax4 = fig.add_subplot(spec[2,1]) + ax4.axis('off') + # 'MSD:' + str(rank_data[0,4]) + ' p-value:' + str(rank_data[0,5])) # --> move to textbox 00 UTC + ##rank_hist_00 =r'$\bf{00 UTC}$\n:' + r'$\bf{MSD}$:' + str(np.round(rank_data[0,4], ndigits)) + '\n' + r'$\bf{ p-value}:' + str(np.round(rank_data[0,5], ndigits)) + rank_hist_00 =r'$\bf{00 UTC}$: ' + r'$\bf{MSD}$: ' + str(np.round(rank_data[0,4], ndigits)) +' '+ r'$\bf{ p-value}$: ' + str(np.round(rank_data[0,5], ndigits)) + + #'MSD:' + str(rank_data[2,4]) + ' p-value:' + str(rank_data[2,5])) --> move to textbox 12 UTC + ##rank_hist_12 =r'$\bf{12 UTC}$\n:' + r'$\bf{MSD}$:' + str(np.round(rank_data[2,4], ndigits)) + '\n' + r'$\bf{ p-value}:' + str(np.round(rank_data[2,5], ndigits)) + rank_hist_12 =r'$\bf{12 UTC}$: ' + r'$\bf{MSD}$: ' + str(np.round(rank_data[2,4], ndigits)) +' '+ r'$\bf{ p-value}$: ' + str(np.round(rank_data[2,5], ndigits)) + rank_hist_text = rank_hist_00 + '\n' + rank_hist_12 + # add text boxes with rank histogram data + ax4.text(0.1,0.4, rank_hist_text, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + + # figure explanation: + expl_data = 'Mean Square Deviation (' + r'$\bf{MSD}$' + ') relative to flat histogram\n' + r'$\bf{p-value}$' + ' of MSD relative to sampling distribution\n' + expl_quantile='Quantiles are estimated based on observations at the station ' + station_id + ' (period: 1978-2020)\n' + expl_model='Model data from experiment (' + exp_id + ') interpolated to location of station ' + station_id + + fig_text = expl_data + expl_quantile+expl_model + + # axis for figure text: + ax5 = fig.add_subplot(spec[3,:]) + ax5.axis('off') + # add text box with explanation to the bottom of the figure + ax5.text(0.1, 0.4,fig_text, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs4) if savefig: plt.savefig(fig_name, dpi=300) -# plt.show( ) - - def read_quantiles_data(quantile_sel_file): ''' Reads quantile data from text file and returns a structured numpy array @@ -165,8 +299,8 @@ def main(config_file): # data to be plotted variable = config['data']['variable'] station_id = config['data']['station_id'] - longitude = config['data']['longitude'] - latitude = config['data']['latitude'] + lon = config['data']['longitude'] + lat = config['data']['latitude'] year_beg = int(config['data']['year_beg']) year_end = int(config['data']['year_end']) month_beg = config['data']['month_beg'] @@ -186,9 +320,7 @@ def main(config_file): rank_hist_data = read_rank_data(rank_hist_file) ## Plot data ## - fig_title = 'Station: ' + station_id + ', ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end + '\nLat=' + latitude + ' Lon=' + longitude - plot_data(quantiles_data, rank_hist_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, True) - + plot_data(quantiles_data, rank_hist_data, time_series_data, station_id, lon, lat, year_beg, year_end, month_beg, month_end, fig_name, True) if __name__=='__main__': if (len(sys.argv) < 2): diff --git a/RADSOUND/STATRS/python/plot_rank_hist_sum_all_stations.py b/RADSOUND/STATRS/python/plot_rank_hist_sum_all_stations.py deleted file mode 100644 index 2648455dd508a423e076d0a39ddcfa76d5106da0..0000000000000000000000000000000000000000 --- a/RADSOUND/STATRS/python/plot_rank_hist_sum_all_stations.py +++ /dev/null @@ -1,94 +0,0 @@ -#!/usr/bin/env python - -''' -Written by Madeleine Ekblom, 2023-09-12 - -Based on Jouni Räisänen's script for plotting rank histogram summart statistics - -Example: -$ python3 plot_rank_hist_sum_all_stations.py rh_summary_file nstat p01 p1 p5 max_freq_p max_freq_q - -''' - -import numpy as np -import matplotlib.pyplot as plt -from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) -import sys - -#def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname): -#def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q): -def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, utc_term): - ''' - Input: p_freq, q_freq, number of stations, p01, p1, p5, max freq p, max freq q, figname - Description: Plots p and q frequencies and saves the figure with figname - ''' - if utc_term == 'UTC00': - title='00 UTC : Rank Histogram Summary Statistics\n Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01:' + str(p1) + '\nFrequency, p<0.05:' + str(p5) - else: - title='12 UTC : Rank Histogram Summary Statistics\n Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01:' + str(p1) + '\nFrequency, p<0.05:' + str(p5) - - # can also be input arguments if the number of bins are not constant - number_of_p_bins=20 - number_of_q_bins=100 - - fig = plt.figure(figsize=(8,12)) - gs = fig.add_gridspec(2, hspace=0.4) - ax1, ax2 = gs.subplots(sharex=False, sharey=False) - - # Plot p values - ax1.set_title('Normalized p-value frequencies') - ax1.set_xlabel('p-value') - ax1.set_ylim([0,max_freq_p]) - ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq) - ax1.xaxis.set_major_locator(MultipleLocator(2)) - ax1.xaxis.set_minor_locator(MultipleLocator(1)) - p_xlabels = ax1.get_xticks() - ax1.set_xticks(p_xlabels, p_xlabels/number_of_p_bins) - ax1.set_xlim([0,number_of_p_bins]) - - # Plot q values - ax2.set_title('Normalized quantile frequencies') - ax2.set_xlabel('Quantile (%)') - ax2.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq) - ax2.xaxis.set_major_locator(MultipleLocator(10)) - ax2.xaxis.set_minor_locator(MultipleLocator(5)) - ax2.set_ylim([0,max_freq_q]) - ax2.set_xlim([0,number_of_q_bins]) - - fig.suptitle(title) - figname='rank_hist_sumstats.png' - plt.savefig(figname, dpi=300) - -def read_data(rh_summary_file): - ''' - Input: path/to/rh_summary_file - Output: p_freq, q_freq - Description: reads in p frequencies and q frequencies - ''' - - # first line contains a header and first column contains hour - p_freq_all = np.loadtxt(rh_summary_file + '_p-freq.txt', skiprows=1) - q_freq_all = np.loadtxt(rh_summary_file + '_q-freq.txt', skiprows=1) - - # p_freq and q_freq contain data 00, 06, 12, 18, 24 (the last is used in these plots) - p_freq = p_freq_all[-1,1:] - q_freq = q_freq_all[-1,1:] - - return p_freq, q_freq - -if __name__=='__main__': - rh_summary_file = sys.argv[1] - number_of_stations = sys.argv[2] - p01 = sys.argv[3] - p1 = sys.argv[4] - p5 = sys.argv[5] - max_freq_p = float(sys.argv[6]) - max_freq_q = float(sys.argv[7]) - utc_term = sys.argv[8] - #figname = sys.argv[8] - - p_freq, q_freq = read_data(rh_summary_file) - #plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname) - #plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q) - plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, utc_term) - diff --git a/RADSOUND/STATRS/python/plot_rank_hist_sum_temp.py b/RADSOUND/STATRS/python/plot_rank_hist_sum_temp.py new file mode 100644 index 0000000000000000000000000000000000000000..480887c71eddc364dea51f86903712464726ec89 --- /dev/null +++ b/RADSOUND/STATRS/python/plot_rank_hist_sum_temp.py @@ -0,0 +1,151 @@ +#!/usr/bin/env python + +''' +Written by Madeleine Ekblom, Oct 2023 +Updates by M.E, 2023-11-29 +- split title into title and textbox with additional information + +Based on Jouni Räisänen's script for plotting rank histogram summary statistics + +Example: +$ python3 plot_rank_hist_sum_temp.py rh_summary_file nstat p01_00 p1_00 p5_00 p01_12 p1_12 p5_12 p_ylim q_ylim figname + +''' + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) +import sys + +def plot_rank_hist(p_freq_00, p_freq_12, q_freq_00, q_freq_12, number_of_stations, p01_00, p1_00, p5_00, p01_12, p1_12, p5_12, p_ylim, q_ylim, figname): + ''' + Input: p_freq 00/12, q_freq 00/12, number of stations, p01, p1, p5, (00/12), p_ylim, q_ylim, figname + Description: Plots p and q frequencies and saves the figure with figname + ''' + # title text + title='RADSOUND: Rank histogram summary statistics\nfor $\Delta$T = 2t - t850' + + # can also be input arguments if the number of p and q bins are not constant + number_of_p_bins=20 + number_of_q_bins=100 + + # set up figure with 2x2 subplots + ##fig = plt.figure(figsize=(9,12)) + fig = plt.figure(figsize=(10,12)) + spec = fig.add_gridspec(3,2, hspace=0.5) + + # plotting details: + fs1 = 20 # fontsize for title + fs2 = 14 # fontsize for labels, subplot titles + fs3 = 10 # fontsize for figure explanations + + # Plot p values, 00 UTC + ax1 = fig.add_subplot(spec[0,0]) + ax1.set_title('00 UTC\nNormalized p-value frequencies', fontsize=fs2) + ax1.set_xlabel('p-value', fontsize=fs2) + ax1.set_ylabel('Relative frequency', fontsize=fs2) + ax1.set_ylim([0,p_ylim]) + ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq_00) # add +0.5 to center the bars + ax1.xaxis.set_major_locator(MultipleLocator(2)) + ax1.xaxis.set_minor_locator(MultipleLocator(1)) + p_xlabels = ax1.get_xticks() + ax1.set_xticks(p_xlabels, p_xlabels/number_of_p_bins) + ax1.set_xlim([0,number_of_p_bins]) + + # Plot p values, 12 UTC + ax2 = fig.add_subplot(spec[0,1]) + ax2.set_title('12 UTC\nNormalized p-value frequencies', fontsize=fs2) + ax2.set_xlabel('p-value', fontsize=fs2) + ax2.set_ylabel('Relative frequency', fontsize=fs2) + ax2.set_ylim([0,p_ylim]) + ax2.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq_12) # add +0.5 to center the bars + ax2.xaxis.set_major_locator(MultipleLocator(2)) + ax2.xaxis.set_minor_locator(MultipleLocator(1)) + p_xlabels = ax2.get_xticks() + ax2.set_xticks(p_xlabels, p_xlabels/number_of_p_bins) + ax2.set_xlim([0,number_of_p_bins]) + + # Plot q values, 00 UTC + ax3 = fig.add_subplot(spec[1,0]) + ax3.set_title('00 UTC\nNormalized quantile frequencies', fontsize=fs2) + ax3.set_xlabel('Quantile (%)', fontsize=fs2) + ax3.set_ylabel('Relative frequency', fontsize=fs2) + ax3.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq_00) # add +0.5 to center the bars + ax3.xaxis.set_major_locator(MultipleLocator(10)) + ax3.xaxis.set_minor_locator(MultipleLocator(5)) + ax3.set_ylim([0,q_ylim]) + ax3.set_xlim([0,number_of_q_bins]) + ax3.axhline(y=1, c='k', ls='--') # horizontal black dashed line at y=0 + + # Plot q values, 12 UTC + ax4 = fig.add_subplot(spec[1,1]) + ax4.set_title('12 UTC\nNormalized quantile frequencies', fontsize=fs2) + ax4.set_xlabel('Quantile (%)', fontsize=fs2) + ax4.set_ylabel('Relative frequency', fontsize=fs2) + ax4.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq_12) # add +0.5 to center the bars + ax4.xaxis.set_major_locator(MultipleLocator(10)) + ax4.xaxis.set_minor_locator(MultipleLocator(5)) + ax4.set_ylim([0,q_ylim]) + ax4.set_xlim([0,number_of_q_bins]) + ax4.axhline(y=1, c='k', ls='--') # horizontal black dashed line at y=0 + + # add text box with additional information: + ax5 = fig.add_subplot(spec[2,:]) + ax5.axis('off') + + # textbox text + textbox = ' Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001, 00 UTC: ' + str(p01_00) + ' Frequency, p<0.001, 12 UTC: ' + str(p01_12) + '\nFrequency, p<0.01, 00 UTC: ' + str(p1_00) + ' Frequency, p<0.01, 12 UTC: ' + str(p1_12) + '\nFrequency, p<0.05, 00 UTC: ' + str(p5_00) + ' Frequency, p<0.05, 12 UTC: ' + str(p5_12) + + ##ax5.text(0.3, 0.6, textbox, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs2) + ax5.text(0.01, 0.6, textbox, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs2) + + # explanation for figure + text_pvalue= 'Normalized p-values: Frequency distribution of p-values at\nall individual stations, normalized to give an expected value of 1.\n' + text_quantile='Normalized quantiles: Frequency distribution of model data\nrelative to quantiles of observations, averaged over all stations\nand divided by the expected value of 0.01.' + + text_expl = text_pvalue + text_quantile + + # add text box with explanation to the bottom of the figure + ax5.text(0.20, -0.05, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + + # add title to figure + fig.suptitle(title, fontsize=fs1) + + plt.savefig(figname) + +def read_data(rh_summary_file): + ''' + Input: path/to/rh_summary_file + Output: p_freq_00, p_freq_12, q_freq_00, q_freq_12 + Description: reads in p frequencies and q frequencies + ''' + #print(rh_summary_file) + + # first line contains a header and first column contains hour + p_freq_all = np.loadtxt(rh_summary_file + '_p-freq.txt', skiprows=1) + q_freq_all = np.loadtxt(rh_summary_file + '_q-freq.txt', skiprows=1) + + # p_freq and q_freq contain data 00, 12, 24 + p_freq_00 = p_freq_all[0,1:] + p_freq_12 = p_freq_all[1,1:] + q_freq_00 = q_freq_all[0,1:] + q_freq_12 = q_freq_all[1,1:] + + return p_freq_00, p_freq_12, q_freq_00, q_freq_12 + +if __name__=='__main__': + rh_summary_file = sys.argv[1] + number_of_stations = sys.argv[2] + p01_00 = sys.argv[3] + p1_00 = sys.argv[4] + p5_00 = sys.argv[5] + p01_12 = sys.argv[6] + p1_12 = sys.argv[7] + p5_12 = sys.argv[8] + p_ylim = float(sys.argv[9]) + q_ylim = float(sys.argv[10]) + figname = sys.argv[11] + + p_freq_00, p_freq_12, q_freq_00, q_freq_12 = read_data(rh_summary_file) + plot_rank_hist(p_freq_00, p_freq_12, q_freq_00, q_freq_12, number_of_stations, p01_00, p1_00, p5_00, p01_12, p1_12, p5_12, p_ylim, q_ylim, figname) + diff --git a/RADSOUND/STATRS/summary_rank_histograms_all_stations.sh b/RADSOUND/STATRS/summary_rank_histograms_all_stations.sh index ba90835e3883937f404d40ff01b9c9612addad8a..21d5ea18a743366a906cb603f9b44afbd4a775bf 100755 --- a/RADSOUND/STATRS/summary_rank_histograms_all_stations.sh +++ b/RADSOUND/STATRS/summary_rank_histograms_all_stations.sh @@ -134,78 +134,93 @@ echo " Calculating the rank histogram summary statistics ..." text_output='${rh_summary_file}_text_output' nbins_p=20 hour1=0,hour2=24,hour_step=12 - grads_out=.true. + grads_out=.false. l_code_in_char=.false. &END EOF + +############################################################## +echo " Converting txt-file to odb-file ... " +echo " ${rh_summary_file}_p-values.txt to ${rh_summary_file}_p-values.odb " + +# Convert '${rh_summary_file}_p-values.txt' to and .odb file +# in order to facilitate extraction of different UTC times for plotting on maps +# +odb import -d ' ' ${rh_summary_file}_p-values.txt ${rh_summary_file}_p-values.odb + ############################################################## echo " Reading the number of stations and all-UTC p-value frequencies ..." # from '${rh_summary_file}_text_output' file +# These will be used as arguments in the plotting script 'plot_rank_hist_sum_temp.py' + +head -1 ${rh_summary_file}_text_output > line1 +read v1 v2 v3 number_of_stations < line1 +echo 'nstat' $number_of_stations +# +head -3 ${rh_summary_file}_text_output | tail -1 > line1 +read v1 p_ylim < line1 +echo 'p ylim' $p_ylim +# +head -4 ${rh_summary_file}_text_output | tail -1 > line1 +read v1 q_ylim < line1 +echo 'q ylim' $q_ylim # -head -1 ${rh_summary_file}_text_output > line_1 -read v1 v2 v3 number_of_stations < line_1 -head -2 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 v2 v3 v4 number_of_p_bins < line_1 -head -3 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 max_freq_p < line_1 -head -4 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 max_freq_q < line_1 -head -5 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 v2 v3 v4 v5 v6 p01_00 < line_1 -head -6 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 v2 v3 v4 v5 v6 p1_00 < line_1 -head -7 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 v2 v3 v4 v5 v6 p5_00 < line_1 -head -8 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 v2 v3 v4 v5 v6 p01_12 < line_1 -head -9 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 v2 v3 v4 v5 v6 p1_12 < line_1 -head -10 ${rh_summary_file}_text_output | tail -1 > line_1 -read v1 v2 v3 v4 v5 v6 p5_12 < line_1 -numbers_of_p_bins_up=$( echo "(${number_of_p_bins} + 0.5)" | bc ) -echo $( echo "(${number_of_p_bins} + 0.5)" | bc ) > line_1 -read number_of_p_bins_up < line_1 -rm line_1 +head -5 ${rh_summary_file}_text_output | tail -1 > line1 +read v1 v2 v3 v4 v5 v6 p01_00 < line1 +echo 'p01 00' $p01_00 + +head -6 ${rh_summary_file}_text_output | tail -1 > line1 +read v1 v2 v3 v4 v5 v6 p1_00 < line1 +echo 'p1 00' $p1_00 + +head -7 ${rh_summary_file}_text_output | tail -1 > line1 +read v1 v2 v3 v4 v5 v6 p5_00 < line1 +echo 'p5 00' $p5_00 + +head -8 ${rh_summary_file}_text_output | tail -1 > line1 +read v1 v2 v3 v4 v5 v6 p01_12 < line1 +echo 'p01 12' $p01_12 + +head -9 ${rh_summary_file}_text_output | tail -1 > line1 +read v1 v2 v3 v4 v5 v6 p1_12 < line1 +echo 'p1 12' $p1_12 + +head -10 ${rh_summary_file}_text_output | tail -1 > line1 +read v1 v2 v3 v4 v5 v6 p5_12 < line1 +echo 'p5 12' $p5_12 + +rm line1 ###################################################################### # -echo " (1) Plotting - P-Values Summary on 2D map ..." +echo " (1) Plotting - P-Values Summary on 2D map for 00 and 12 UTCs ..." # ###################################################################### -python_script=python/plot_p_values_map.py +odb sql 'select station,longitude,latitude,p_value where hour=00' -i ${rh_summary_file}_p-values.odb -o p_values_00.txt +odb sql 'select station,longitude,latitude,p_value where hour=12' -i ${rh_summary_file}_p-values.odb -o p_values_12.txt + +python_script=python/plot_p_values_map_rank_hist_temp.py -python3 ${python_script} ${rh_summary_file}_p-values.txt ${year1}${month1}-${year2}${month2} -mv p_values.png ${figure_dir}/p-value_map_${variable}_${year1}${month1}-${year2}${month2}.png -echo " Output file: ${figure_dir}/p-value_map_${variable}_${year1}${month1}-${year2}${month2}.png" +figname=${figure_dir}/p-values_rank-histograms_${variable}_${year1}${month1}-${year2}${month2}.png +echo " Output file: $figname" + +python3 ${python_script} p_values_00.txt p_values_12.txt ${year1}${month1}-${year2}${month2} ${figname} + +rm p_values_*.txt ##################################################################### # echo " (2) Plotting - Rank Histogram Summary Statistics for 00 and 12 UTCs ..." # ##################################################################### -utc_term="UTC00" -echo " for 00 UTC - variable : ${variable}" -python_script=python/plot_rank_hist_sum_all_stations.py -echo " number_of_stations : ${number_of_stations}" -echo " p01_00, p1_00, p5_00, max_freq_p, max_freq_q : ${p01_00} ${p1_00} ${p5_00} ${max_freq_p} ${max_freq_q}" - -python3 ${python_script} ${rh_summary_file} ${number_of_stations} ${p01_00} ${p1_00} ${p5_00} ${max_freq_p} ${max_freq_q} ${utc_term} ${figname} -mv rank_hist_sumstats.png ${figure_dir}/rank-hist-00utc-sumstats_${variable}_${year1}${month1}-${year2}${month2}.png -echo " Output file for 00 UTC: ${figure_dir}/rank-hist-00utc-sumstats_${variable}_${year1}${month1}-${year2}${month2}.png" - -##################################################################### -utc_term="UTC12" -echo " for 12 UTC - variable : ${variable}" -python_script=python/plot_rank_hist_sum_all_stations.py -echo " number_of_stations : ${number_of_stations}" -echo " p01_12, p1_12, p5_12, max_freq_p, max_freq_q : ${p01_12} ${p1_12} ${p5_12} ${max_freq_p} ${max_freq_q}" +python_script=python/plot_rank_hist_sum_temp.py -python3 ${python_script} ${rh_summary_file} ${number_of_stations} ${p01_12} ${p1_12} ${p5_12} ${max_freq_p} ${max_freq_q} ${utc_term} ${figname} +figname=${figure_dir}/summary_plot_rank-histograms_${variable}_${year1}${month1}-${year2}${month2}.png +echo " Output file: $figname" -mv rank_hist_sumstats.png ${figure_dir}/rank-hist-12utc-sumstats_${variable}_${year1}${month1}-${year2}${month2}.png -echo " Output file for 12 UTC: ${figure_dir}/rank-hist-12utc-sumstats_${variable}_${year1}${month1}-${year2}${month2}.png" +python3 ${python_script} ${rh_summary_file} ${number_of_stations} ${p01_00} ${p1_00} ${p5_00} ${p01_12} ${p1_12} ${p5_12} ${p_ylim} ${q_ylim} ${figname} ###################################################################### # Delete files that are not needed any more diff --git a/RADSOUND/main_radsound.sh b/RADSOUND/main_radsound.sh index 31234ea3cca92aac14ab3a1649fd358e9c078510..aa883cbba4fe4e7b23eaaded1585a50fd0f0a9e9 100755 --- a/RADSOUND/main_radsound.sh +++ b/RADSOUND/main_radsound.sh @@ -97,10 +97,8 @@ echo " STEP #1 : PRE-PROCESSING MOD DATA EXTRACTED WITH GSV_INTERFACE ..." echo " SCRIPT --- gsv_radsound_mod_data.sh" echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" -####### ./gsv_radsound_mod_data.sh - echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" echo " ..... STEP 1 - COMPLETED ....." echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" diff --git a/README.md b/README.md index 44fdc1c0a852e17952369656d6cad1180ad9de57..e33e6e6ca1b05292db97a0cc85829b69a9bf0a80 100644 --- a/README.md +++ b/README.md @@ -212,15 +212,23 @@ sys.exit(0) ``` -**Video-recording of installaion and demo-run:** ... (coming soon) ... +# Video-recording & slides of installaion and demo-run: ... (coming soon) ... + +
+ OBSALL Apps installation +
OBSALL Apps installation
+
+ + -
- -``` -... - -``` -
# Disclaimer The OBSALL package is in a developement phase by the University of Helsinki team, led by Heikki Järvinen (heikki.j.jarvinen@helsinki.fi). The team includes also Jouni Räisänen (jouni.raisanen@helsinki.fi), Lauri Tuppi (lauri.tuppi@helsinki.fi), Madeleine Ekblom (madeleine.ekblom@helsinki.fi), and Alexander Mahura (alexander.mahura@helsinki.fi). Some features are still not implemented and you may expect bugs. For feedback and issue reporting, feel free to open an issue in: https://earth.bsc.es/gitlab/digital-twins/de_340/obsall/-/issues diff --git a/SYNOP/STATS/python/plot_p_values_map.py b/SYNOP/STATS/python/plot_p_values_map.py index 7b6e3527bb32cf416e20330327ff36769ca5f2bb..4c17da7b5b1e6d0e9ad8e494230425161e861f2d 100644 --- a/SYNOP/STATS/python/plot_p_values_map.py +++ b/SYNOP/STATS/python/plot_p_values_map.py @@ -2,6 +2,8 @@ ''' Written by Madeleine Ekblom, 2023-04-06 +updates 2023-11-27 +- Option of moving legend into plot when plotting on a global map(lower left corner) Based on Jouni Räisänen's script for plotting p values on the map of Finland @@ -38,43 +40,95 @@ def plot_p_values(p_values): Input: numpy array Output: saves a figure in the running directory ''' - - lon_min=19 - lon_max=32 - lat_min=59 - lat_max=71 - + # add more variables and units when needed + variable_shortnames={'39': '2t'} + variable_units={'39': '$^\circ$C'} + + # variables to add, shortnames and units need to be checked + # total amount of clouds, 91 + # sea level pressure, 108 + # 1-hour precipitation, 80 + # relative humidity, 58 + # 10-minute precipitation intensity, 999 + # snow depth, 71 + # 2m temp, 39 --> 2t, degrees Celcius + # dew point temperature, 40 + # visibility, 6 + # wind direction, 111 + # max wind gust in last 10 mins,261 + # surface wind speed, 221 + + variable='39' # change to be an input value (as string) + + # if data only from SYNOP stations in Finland: + only_Finland_data=True + + # Define area of map: + if only_Finland_data: + # Finland: + lon_min=19 + lon_max=32 + lat_min=59 + lat_max=71 + else: + # Global: + lon_min=0 + lon_max=359 + lat_min=-90 + lat_max=90 + + # define colours and limits for the colours: colors = ['darkviolet', 'blue', 'skyblue', 'lime', 'yellow', 'orange', 'red', 'grey'] limits = [0,0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,1] - fig = plt.figure(figsize=(12,8)) + # create figure and add a subplot with a map: + fig = plt.figure(figsize=(8,12)) ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree()) ax.coastlines() ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree()) ax.gridlines(draw_labels=True) - # basic plot, not changing colors of p: - #ax.scatter(p_values['lon'], p_values['lat'], c=p_values['p']) - + # Plotting details: + fs1=20 # fontsize for title + fs2=14 # fontsize for legend + fs3=10 # fontisze for text box explanation + #loop over limits: legend_elements = [] for n in range(8): + # find indices p_ind = ((p_values['p'] > limits[n]) & (p_values['p'] <= limits[n+1])) + # scatter plot ax.scatter(p_values['lon'][p_ind], p_values['lat'][p_ind], c=colors[n]) + # add to legend_elements, color and values legend_elements.append(Line2D([0], [0], marker='o', color='w',markerfacecolor=colors[n], label=str(limits[n])+'-'+str(limits[n+1]))) - # legend table with customized values - - #print(legend_elements) - - ax.legend(handles=legend_elements, loc='center left', bbox_to_anchor=(1.1, 0.5)) - ax.set_title('P-values for T2m quantiles in ' + timestring, fontsize=20) - + # explanation for figure: + text_expl ='p-value: the probability that the MSD statistics for the quantile rank\nhistogram would be exceeded by pure chance\nwhen comparing a perfect model simulation with observations.' + + # add legend table with customized values + if only_Finland_data: + # If Finland map, put legend outside map + ax.legend(handles=legend_elements, loc='center left', fontsize=fs3) #, bbox_to_anchor=(1.1, 0.5), fontsize=fs2) + # add text box with explanation to the bottom of the figure: 0.1, -0.25 + ax.text(0.1, -0.25, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3, transform=ax.transAxes) + else: + # If global map, put Legend in the lower left corner + ax.legend(handles=legend_elements, loc='lower left', fontsize=fs3) + # add text box with explanation to the bottom of the figure, 0.1, -0.45 + ax.text(0.1, -0.45, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3, transform=ax.transAxes) + + # set title: + fig.suptitle('SYNOP stations: p-values for ' + variable_shortnames[variable] +'\nquantiles in ' + timestring, fontsize=fs1) + + # save figure plt.savefig('p_values.png', dpi=300) -# plt.show() def main(p_values_file): + ''' + Main file: (1) reads in p values from file, (2) plot p values on a map + ''' p_values_data = read_data(p_values_file) diff --git a/SYNOP/STATS/python/plot_quantiles_rankhist.py b/SYNOP/STATS/python/plot_quantiles_rankhist.py index b5fd9958ffa4fe85ca38d8a56ea5c8962e8ed7ef..aec03020e2c76c48bce952812d05a1fd0ac46ac8 100644 --- a/SYNOP/STATS/python/plot_quantiles_rankhist.py +++ b/SYNOP/STATS/python/plot_quantiles_rankhist.py @@ -2,6 +2,8 @@ ''' Written by Madeleine Ekblom, 2023-05-10 +updates 2023-11-27: +- visual changes of the plot: labels, legends , title, and colours updated Based on Jouni Räisänen's script for plotting quantiles/time series data (00, 06, 12, 18 UTC) and rank histogram (24UTC=combination of 00,06,12,and 18 UTC data) @@ -35,35 +37,109 @@ import sys import configparser import datetime from matplotlib.dates import drange, MonthLocator, DateFormatter +import matplotlib.lines as mlines +import matplotlib.patches as mpatches +import matplotlib.ticker as ticker config = configparser.ConfigParser() -#plot_data(quantiles_data, rank_hist_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, True) -def plot_data(quantiles_data, rank_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, savefig=False): +def plot_data(quantiles_data, rank_data, time_series_data, variable, station_id, lon, lat, year_beg, year_end, month_beg, month_end, figname, savefig=False): ''' - quantiles data give the lines and the time series data give the dots + Plot data: quantiles data give the lines and the time series data give the dots ''' + # add more variables and units when needed + variable_shortnames={'39': '2t'} + variable_units={'39': '$^\circ$C'} + + # variables to add, shortnames and units need to be checked + # total amount of clouds, 91 + # sea level pressure, 108 + # 1-hour precipitation, 80 + # relative humidity, 58 + # 10-minute precipitation intensity, 999 + # snow depth, 71 + # 2m temp, 39 --> 2t, degrees Celcius + # dew point temperature, 40 + # visibility, 6 + # wind direction, 111 + # max wind gust in last 10 mins,261 + # surface wind speed, 221 + + # experiment id: now use 'test data' - to be replaced with real exp_id (as input argument) + exp_id = 'test data' + # set up date array consisting of all days in beg year: days = np.arange(1,366) -# day1 = datetime.date(year_beg, 1, 1) -# day2 = datetime.date(year_beg+1, 1, 1) + #day1 = datetime.date(year_beg, 1, 1) + #day2 = datetime.date(year_beg+1, 1, 1) day1 = datetime.date(2001, 1, 1) day2 = datetime.date(2002, 1, 1) delta = datetime.timedelta(days=1) dates = drange(day1, day2, delta) - + # arrays for quantile names, axes titles, hours to plot q_names = ['q01', 'q10', 'q25', 'q50', 'q75', 'q90', 'q99'] sub_titles = ['00 UTC', '06 UTC', '12 UTC', '18 UTC'] hours = [0, 6, 12, 18] - + + # plotting details: + alpha = [0.5, 0.3, 0.1] # transparency for shaded area plotting + color='blue' # color for shading + linecolor='red' # color for solid line + scatter_size=5 # size for scattered dots + fs1=20 # fontsize for title + fs2=14 # fontsize for labels, subplot titles + fs3 = 12 # fontsize for + fs4=10 # fontsize for extra text box, legend + ndigits = 3 # round off to ndigits (e.g., lon, lat) + + # use variable to get variable and corresponding unit, used as ylabel in plots: + ylabel = variable_shortnames[variable] + ' (' + variable_units[variable] + ')' + + # legend elemnents for quantile data, one legend for model and one legend for observations + legend_elements1 = [mlines.Line2D([],[], marker='.', linestyle='', color='k', label='Model data (' + exp_id + ')')] + legend_elements2= [mpatches.Patch(color=color, alpha=alpha[0], label='25-50%'), + mpatches.Patch(color=color, alpha=alpha[1], label='10-25%'), + mpatches.Patch(color=color, alpha=alpha[2], label='1-10%'), + mlines.Line2D([],[], color=linecolor, label='50%'), + mpatches.Patch(color=color, alpha=alpha[0], label='50-75%'), + mpatches.Patch(color=color, alpha=alpha[1], label='75-90%'), + mpatches.Patch(color=color, alpha=alpha[2], label='90-99%')] + + # round lon, lat to three digits: + lon = np.round(float(lon), ndigits) + lat = np.round(float(lat), ndigits) + + # define lon/lat text strings for title of the form lat: XX N/S, lat: XX E/W + if (lat > 0): + # if lat positive -> N + lat_text = ' (lat: ' + str(lat) + '$^\circ$N, ' + else: + # if negative -> S + lat_text = ' (lat ' + str(-lat) + '$^\circ$S, ' + + if (0 <= lon) and (lon <= 180): + # lon btw 0 an 180 -> E + lon_text = ' lon: ' + str(lon) + '$^\circ$E)' + else: + if lon < 0: + # lon < 0 --> - lon W + lon_text = ' lon: ' + str(-lon) + '$^\circ$W)' + else: + # 180 < lon < 360 --> - (lon - 360) W + lon = -(lon - 360) + lon_text = ' lon: ' + str(lon) + '$^\circ$W)' + + # define figure title: + fig_title = 'SYNOP station: ' + station_id + lat_text + lon_text + '\nModel data from ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end + # calculate number of years: nyears = year_end - year_beg + 1 # set up figure fig = plt.figure(figsize=(10,10), layout='constrained') - spec = fig.add_gridspec(3,2) - fig.suptitle(fig_title, fontsize=20) + spec = fig.add_gridspec(4,2) + fig.suptitle(fig_title, fontsize=fs1) # counter c = 0 @@ -71,40 +147,95 @@ def plot_data(quantiles_data, rank_data, time_series_data, fig_title, fig_name, # plot quantiles/time series for times 00, 06, 12, 18: for i in range(2): for j in range(2): - # set up axis: title, xaxis + # set up axis: title, xaxis, x- and y-labels ax = fig.add_subplot(spec[i,j]) - ax.set_title(sub_titles[c]) + ax.set_title(sub_titles[c], fontsize=fs2) ax.set_xlim(dates[0], dates[-1]) ax.xaxis.set_major_locator(MonthLocator()) ax.xaxis.set_major_formatter(DateFormatter('%b')) - + ax.set_xlabel('Time', fontsize=fs2) + ax.set_ylabel(ylabel, fontsize=fs2) + # quantile data: # find quantile data hour = 0, 6, 12, 18 qh_ind = (quantiles_data['hour'][:] == hours[c]) - for q in q_names: - q_data = quantiles_data[q][qh_ind] - ax.plot(dates, q_data[:365], 'k-') # --> change q only contains one year of data!! + + # get quantile data, qXX_data = quantile XX data at hour qh_ind + q01_data = quantiles_data[q_names[0]][qh_ind] + q10_data = quantiles_data[q_names[1]][qh_ind] + q25_data = quantiles_data[q_names[2]][qh_ind] + q50_data = quantiles_data[q_names[3]][qh_ind] + q75_data = quantiles_data[q_names[4]][qh_ind] + q90_data = quantiles_data[q_names[5]][qh_ind] + q99_data = quantiles_data[q_names[6]][qh_ind] + + # solid line at q50: + ax.plot(dates, q50_data, c=linecolor, ls='-') + + # shading will be the most opaque close to 50 and more transparent close to 1/99 + + # fill area btw q25-q50, q50-q75: alpha = alpha[0] + ax.fill_between(dates, q25_data[:], q50_data[:], color=color, alpha=alpha[0]) + ax.fill_between(dates, q50_data[:], q75_data[:], color=color, alpha=alpha[0]) + + # fill area btw q10-q25, q75-q90, alpha=alpha[1] + ax.fill_between(dates, q10_data[:], q25_data[:], color=color, alpha=alpha[1]) + ax.fill_between(dates, q75_data[:], q90_data[:], color=color, alpha=alpha[1]) + + # fill area btw q01-q10, q90-q99, alpha=alpha[2] + ax.fill_between(dates, q01_data[:], q10_data[:], color=color, alpha=alpha[2]) + ax.fill_between(dates, q90_data[:], q99_data[:], color=color, alpha=alpha[2]) + # plot time series data: # find time series where hour = 0, 6, 12, 18 th_ind = (time_series_data['hour'][:] == hours[c]) t_data = time_series_data['value'][th_ind] # all years + for n in range(nyears): - ax.scatter(dates, t_data[n*365:(n+1)*365], marker='.') + # plot data with black, small dots + ax.scatter(dates, t_data[n*365:(n+1)*365], s=scatter_size, c='k', marker='.') c = c + 1 # plot rank histogram data: - ax2 = fig.add_subplot(spec[2,:]) - ax2.set_title('Rank histogram\n' + 'MSD:' + str(rank_data[4,4]) + ' p-value:' + str(rank_data[4,5])) - ax2.bar(np.arange(100), rank_data[4,6:]*100) + ax2 = fig.add_subplot(spec[2,0]) + ax2.set_title('Rank histogram', fontsize=fs2) + ax2.bar(np.arange(100)+0.5, rank_data[4,6:]*100) # add + 0.5 to center bars ax2.set_xlim(0,100) - ax2.axhline(y=1) + ax2.set_xlabel('Percentile', fontsize=fs2) + ax2.set_ylabel('Relative frequency ', fontsize=fs2) + ax2.axhline(y=1, c='k', ls='--') # horizontal line at y=1.0 + ax2.xaxis.set_major_locator(ticker.FixedLocator([10,25,50,75,90])) + + # add a sixth axis for legends and text data: + ax3 = fig.add_subplot(spec[2,1]) + ax3.axis('off') + + # legend with two columns at the upper centre of the figure: + leg1 = ax3.legend(loc='upper center', handles=legend_elements1, fontsize=fs4) + leg2 = ax3.legend(loc='lower center', handles=legend_elements2, ncols=2, fontsize=fs4, title='Quantiles for observations', title_fontsize=fs4) + ax3.add_artist(leg1) + ax3.add_artist(leg2) + + # add text boxes with (1) rank histogram data and (2) explanation + rank_hist_text = r'$\bf{MSD}$:' + str(np.round(rank_data[4,4], ndigits)) + '\n' + r'$\bf{ p-value}$:' + str(np.round(rank_data[4,5], ndigits)) + + # add text boxes to axis in right lower corner + ax3.text(0.35, -0.3, rank_hist_text, bbox={'pad': 10, 'facecolor':'white'}, fontsize=fs3) + # figure explanation: explaining what quantiles and model data are + expl_msd_pvalue = 'Mean Square Deviation (' + r'$\bf{MSD}$' + ' relative to flat histogram\n' + r'$\bf{p-value}$' + ' of MSD relative to sampling distribution\n' + expl_quantile='Quantiles are estimated based on observations at the station ' + station_id +' (period: 2010-2021)\n' + expl_model='Model data from experiment (' + exp_id + ') interpolated to location of this station ' + fig_text = expl_msd_pvalue + expl_quantile + expl_model + + # add subplot at the bottom of the figure, and add figure explanation as text to this subplot + ax4 = fig.add_subplot(spec[3,:]) + ax4.axis('off') + ax4.text(0.1, 0.4, fig_text, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs4) + + # save figure: if savefig: - plt.savefig(fig_name, dpi=300) - - #plt.savefig(fig_name, dpi=300) - -# plt.show( ) + plt.savefig(figname, dpi=300) def read_quantiles_data(quantile_sel_file): @@ -161,8 +292,8 @@ def main(config_file): # data to be plotted variable = config['data']['variable'] station_id = config['data']['station_id'] - longitude = config['data']['longitude'] - latitude = config['data']['latitude'] + lon = config['data']['longitude'] + lat = config['data']['latitude'] year_beg = int(config['data']['year_beg']) year_end = int(config['data']['year_end']) month_beg = config['data']['month_beg'] @@ -174,7 +305,7 @@ def main(config_file): rank_hist_file = config['input']['rank_hists'] # output files - fig_name = config['output']['fig_name'] + figname = config['output']['fig_name'] ## Read input data ## time_series_data = read_time_series_data(sim_file) @@ -182,8 +313,7 @@ def main(config_file): rank_hist_data = read_rank_data(rank_hist_file) ## Plot data ## - fig_title = 'Station: ' + station_id + ', ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end + '\nLat=' + latitude + ' Lon=' + longitude - plot_data(quantiles_data, rank_hist_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, True) + plot_data(quantiles_data, rank_hist_data, time_series_data, variable, station_id, lon, lat, year_beg, year_end, month_beg, month_end, figname, True) if __name__=='__main__': diff --git a/SYNOP/STATS/python/plot_rank_hist_sum_all_stations.py b/SYNOP/STATS/python/plot_rank_hist_sum_all_stations.py index 602caac02c2942ceacec6554b58eac26fb107971..10ebd3235cdfecdd53cf1eebf1339d3a8f5d9b5e 100644 --- a/SYNOP/STATS/python/plot_rank_hist_sum_all_stations.py +++ b/SYNOP/STATS/python/plot_rank_hist_sum_all_stations.py @@ -2,6 +2,8 @@ ''' Written by Madeleine Ekblom, 2023-09-12 +updates 2023-11-27: +- moved text from title to a new text box Based on Jouni Räisänen's script for plotting rank histogram summart statistics @@ -15,27 +17,57 @@ import matplotlib.pyplot as plt from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) import sys -#def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname): def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q): ''' - Input: p_freq, q_freq, number of stations, p01, p1, p5, max freq p, max freq q, figname + Input: p_freq, q_freq, number of stations, p01, p1, p5, max freq p, max freq q Description: Plots p and q frequencies and saves the figure with figname ''' - title='Rank Histogram Summary Statistics\n Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01:' + str(p1) + '\nFrequency, p<0.05:' + str(p5) + # add more variables and units when needed + variable_shortnames={'39': '2t'} + variable_units={'39': '$^\circ$C'} - # can also be input arguments if the number of bins are not constant + # variables to add, shortnames and units need to be checked + # total amount of clouds, 91 + # sea level pressure, 108 + # 1-hour precipitation, 80 + # relative humidity, 58 + # 10-minute precipitation intensity, 999 + # snow depth, 71 + # 2m temp, 39 --> 2t, degrees Celcius + # dew point temperature, 40 + # visibility, 6 + # wind direction, 111 + # max wind gust in last 10 mins,261 + # surface wind speed, 221 + + variable='39' # change to be an input value, as string + + # title + title='SYNOP: Rank histogram summary statistics for ' + variable_shortnames[variable] + + # text box information + text_box = 'Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01: ' + str(p1) + '\nFrequency, p<0.05: ' + str(p5) + + # these can also be input arguments if the number of bins are not constant number_of_p_bins=20 number_of_q_bins=100 + # set up figure fig = plt.figure(figsize=(8,12)) - gs = fig.add_gridspec(2, hspace=0.4) - ax1, ax2 = gs.subplots(sharex=False, sharey=False) + gs = fig.add_gridspec(3, hspace=0.5) + ax1, ax2, ax3 = gs.subplots(sharex=False, sharey=False) + + # plotting details + fs1=20 # fontsize for title + fs2=14 # fontsize for labels, subplot titles, textbox + fs3=10 # fontsize for text explanation # Plot p values - ax1.set_title('Normalized p-value frequencies') - ax1.set_xlabel('p-value') + ax1.set_title('Normalized p-value frequencies', fontsize=fs2) + ax1.set_xlabel('p-value', fontsize=fs2) + ax1.set_ylabel('Relative frequency', fontsize=fs2) ax1.set_ylim([0,max_freq_p]) - ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq) + ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq) # add +0.5 to center the bars ax1.xaxis.set_major_locator(MultipleLocator(2)) ax1.xaxis.set_minor_locator(MultipleLocator(1)) p_xlabels = ax1.get_xticks() @@ -43,15 +75,34 @@ def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, ax1.set_xlim([0,number_of_p_bins]) # Plot q values - ax2.set_title('Normalized quantile frequencies') - ax2.set_xlabel('Quantile (%)') - ax2.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq) + ax2.set_title('Normalized quantile frequencies', fontsize=fs2) + ax2.set_xlabel('Quantile (%)', fontsize=fs2) + ax2.set_ylabel('Relative frequency', fontsize=fs2) + ax2.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq) # add +0.5 to center the bars ax2.xaxis.set_major_locator(MultipleLocator(10)) ax2.xaxis.set_minor_locator(MultipleLocator(5)) ax2.set_ylim([0,max_freq_q]) ax2.set_xlim([0,number_of_q_bins]) - fig.suptitle(title) + # remove axis from ax3 + ax3.axis('off') + + # add text box with additional infromation: + ax3.text(0.3, 0.6, text_box, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs2) + + # explanation for figure + text_pvalue= 'Normalized p-values: Frequency distribution of p-values at\nall individual stations, normalized to give an expected value of 1.\n' + text_quantile='Normalized quantiles: Frequency distribution of model data\nrelative to quantiles of observations, averaged over all stations\nand divided by the expected value of 0.01.' + + text_expl = text_pvalue + text_quantile + + # add text box with explanation to the bottom of the figure + ax3.text(0.1, -0.05, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + + # add title to figure + fig.suptitle(title, fontsize=fs1) + + # save figure figname='rank_hist_sumstats.png' plt.savefig(figname, dpi=300) @@ -80,9 +131,7 @@ if __name__=='__main__': p5 = sys.argv[5] max_freq_p = float(sys.argv[6]) max_freq_q = float(sys.argv[7]) - #figname = sys.argv[8] p_freq, q_freq = read_data(rh_summary_file) - #plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname) plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q) diff --git a/SYNOP/STATS/python/plot_rank_hist_sum_all_stations_v2.py b/SYNOP/STATS/python/plot_rank_hist_sum_all_stations_v2.py deleted file mode 100644 index e31630ae90ef0bc4fa5b4ea28a76d09f8bbca899..0000000000000000000000000000000000000000 --- a/SYNOP/STATS/python/plot_rank_hist_sum_all_stations_v2.py +++ /dev/null @@ -1,85 +0,0 @@ -#!/usr/bin/env python - -''' -Written by Madeleine Ekblom, 2023-09-12 - -Based on Jouni Räisänen's script for plotting rank histogram summart statistics - -Example: -$ python3 plot_rank_hist_sum_all_stations.py rh_summary_file nstat p01 p1 p5 max_freq_p max_freq_q - -''' - -import numpy as np -import matplotlib.pyplot as plt -from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) -import sys - -def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname): - ''' - Input: p_freq, q_freq, number of stations, p01, p1, p5, max freq p, max freq q, figname - Description: Plots p and q frequencies and saves the figure with figname - ''' - title='Rank histogram summary statistics\n Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01: ' + str(p1) + '\nFrequency, p<0.05: ' + str(p5) - - # can also be input arguments if the number of bins are not constant - number_of_p_bins=20 - number_of_q_bins=100 - - fig = plt.figure(figsize=(8,12)) - gs = fig.add_gridspec(2, hspace=0.4) - ax1, ax2 = gs.subplots(sharex=False, sharey=False) - - # Plot p values - ax1.set_title('Normalized p-value frequencies') - ax1.set_xlabel('p-value') - ax1.set_ylim([0,max_freq_p]) - ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq) - ax1.xaxis.set_major_locator(MultipleLocator(2)) - ax1.xaxis.set_minor_locator(MultipleLocator(1)) - p_xlabels = ax1.get_xticks() - ax1.set_xticks(p_xlabels, p_xlabels/number_of_p_bins) - ax1.set_xlim([0,number_of_p_bins]) - - # Plot q values - ax2.set_title('Normalized quantile frequencies') - ax2.set_xlabel('Quantile (%)') - ax2.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq) - ax2.xaxis.set_major_locator(MultipleLocator(10)) - ax2.xaxis.set_minor_locator(MultipleLocator(5)) - ax2.set_ylim([0,max_freq_q]) - ax2.set_xlim([0,number_of_q_bins]) - - fig.suptitle(title) - plt.savefig(figname) - -def read_data(rh_summary_file): - ''' - Input: path/to/rh_summary_file - Output: p_freq, q_freq - Description: reads in p frequencies and q frequencies - ''' - - # first line contains a header and first column contains hour - p_freq_all = np.loadtxt(rh_summary_file + '_p-freq.txt', skiprows=1) - q_freq_all = np.loadtxt(rh_summary_file + '_q-freq.txt', skiprows=1) - - # p_freq and q_freq contain data 00, 06, 12, 18, 24 (the last is used in these plots) - p_freq = p_freq_all[-1,1:] - q_freq = q_freq_all[-1,1:] - - return p_freq, q_freq - -if __name__=='__main__': - rh_summary_file = sys.argv[1] - number_of_stations = sys.argv[2] - p01 = sys.argv[3] - p1 = sys.argv[4] - p5 = sys.argv[5] - max_freq_p = float(sys.argv[6]) - max_freq_q = float(sys.argv[7]) - figname = sys.argv[8] - - p_freq, q_freq = read_data(rh_summary_file) - plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname) -