From f3266dc1d28aec5c7f5430388800a3db7f0a8a74 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Fri, 11 Oct 2024 16:55:17 -0500 Subject: [PATCH 01/10] ADD: New example showing how to use ACT to work with ARM and Satellite data --- act/plotting/xsectiondisplay.py | 46 +++++++++++++++++-------- examples/plotting/plot_satellite.py | 52 +++++++++++++++++++++++++++++ 2 files changed, 84 insertions(+), 14 deletions(-) create mode 100644 examples/plotting/plot_satellite.py diff --git a/act/plotting/xsectiondisplay.py b/act/plotting/xsectiondisplay.py index 0538e8d08d..cc320e28cd 100644 --- a/act/plotting/xsectiondisplay.py +++ b/act/plotting/xsectiondisplay.py @@ -146,13 +146,14 @@ def set_yrng(self, yrng, subplot_index=(0,)): def plot_xsection( self, - dsname, - varname, + field, + dsname=None, x=None, y=None, subplot_index=(0,), sel_kwargs=None, isel_kwargs=None, + set_title=None, **kwargs, ): """ @@ -162,11 +163,10 @@ def plot_xsection( Parameters ---------- - dsname : str or None - The name of the datastream to plot from. Set to None to have - ACT attempt to automatically detect this. - varname : str + field : str The name of the variable to plot. + dsname : str or None + The name of the datastream to plot from. x : str or None The name of the x coordinate variable. y : str or None @@ -181,6 +181,8 @@ def plot_xsection( to slice datasets, see the documentation on :func:`xarray.DataArray.sel`. isel_kwargs : dict The keyword arguments to pass into :py:func:`xarray.DataArray.sel` + set_title : str + Title for the plot **kwargs : keyword arguments Additional keyword arguments will be passed into :func:`xarray.DataArray.plot`. @@ -220,9 +222,9 @@ def plot_xsection( y_coord_dim = temp_ds[y].dims[0] coord_list[y] = y_coord_dim new_ds = data_utils.assign_coordinates(temp_ds, coord_list) - my_dataarray = new_ds[varname] + my_dataarray = new_ds[field] else: - my_dataarray = temp_ds[varname] + my_dataarray = temp_ds[field] coord_keys = [key for key in my_dataarray.coords.keys()] # X-array will sometimes shorten latitude and longitude variables @@ -255,28 +257,42 @@ def plot_xsection( self.set_xrng(xrng, subplot_index) yrng = self.axes[subplot_index].get_ylim() self.set_yrng(yrng, subplot_index) + + if set_title is None: + if 'long_name' in self._ds[dsname][field].attrs: + set_title = self._ds[dsname][field].attrs['long_name'] + plt.title(set_title) + del temp_ds return self.axes[subplot_index] def plot_xsection_map( - self, dsname, varname, subplot_index=(0,), coastlines=True, background=False, **kwargs + self, + field, + dsname=None, + subplot_index=(0,), + coastlines=True, + background=False, + set_title=None, + **kwargs, ): """ Plots a cross section of 2D data on a geographical map. Parameters ---------- - dsname : str or None - The name of the datastream to plot from. Set to None - to have ACT attempt to automatically detect this. - varname : str + field : str The name of the variable to plot. + dsname : str or None + The name of the datastream to plot from. subplot_index : tuple The index of the subplot to plot inside. coastlines : bool Set to True to plot the coastlines. background : bool Set to True to plot a stock image background. + set_title : str + Title for the plot **kwargs : keyword arguments Additional keyword arguments will be passed into :func:`act.plotting.XSectionDisplay.plot_xsection` @@ -292,7 +308,9 @@ def plot_xsection_map( 'Cartopy needs to be installed in order to plot ' + 'cross sections on maps!' ) self.set_subplot_to_map(subplot_index) - self.plot_xsection(dsname, varname, subplot_index=subplot_index, **kwargs) + self.plot_xsection( + field, dsname=dsname, subplot_index=subplot_index, set_title=set_title, **kwargs + ) xlims = self.xrng[subplot_index].flatten() ylims = self.yrng[subplot_index].flatten() self.axes[subplot_index].set_xticks(np.linspace(round(xlims[0], 0), round(xlims[1], 0), 10)) diff --git a/examples/plotting/plot_satellite.py b/examples/plotting/plot_satellite.py new file mode 100644 index 0000000000..8303b594a9 --- /dev/null +++ b/examples/plotting/plot_satellite.py @@ -0,0 +1,52 @@ +""" +Using ACT for Satellite data +---------------------------- + +Simple example for working with satellite data. +It is recommended that users explore other +satellite specific libraries suchas SatPy. + +Author: Adam Theisen +""" + +import act +import matplotlib.pyplot as plt +import numpy as np +from arm_test_data import DATASETS + +# Read in VISST Data +files = DATASETS.fetch('enavisstgridm11minnisX1.c1.20230307.000000.cdf') +ds = act.io.read_arm_netcdf(files) + +# Plot up the VISST cloud percentage using XSectionDisplay +display = act.plotting.XSectionDisplay(ds, figsize=(10, 8)) +display.plot_xsection_map( + 'cloud_percentage', x='longitude', y='latitude', isel_kwargs={'time': 0, 'cld_type': 0} +) +plt.show() + +# Download ARM TSI Data +files = DATASETS.fetch('enatsiskycoverC1.b1.20230307.082100.cdf') +ds_tsi = act.io.read_arm_netcdf(files) +ds_tsi = ds_tsi.where(ds_tsi.percent_opaque > 0) + +# Set coordinates to extra data for ENA +ena_lat = 39.091600 +ena_lon = 28.025700 + +lat = ds['lat'].values +lon = ds['lon'].values + +# Find the nearest pixel for the satellite data and extract it +lat_ind = np.argmin(np.abs(lat - ena_lat)) +lon_ind = np.argmin(np.abs(lon - ena_lon)) + +ds_new = ds.isel(lat=lat_ind, lon=lon_ind, cld_type=0) + +# Plot the comparison using TimeSeriesDisplay +display = act.plotting.TimeSeriesDisplay({'Satellite': ds_new, 'ARM': ds_tsi}, figsize=(15, 8)) +display.plot('cloud_percentage', dsname='Satellite', label='VISST Cloud Percentage') +display.plot('percent_opaque', dsname='ARM', label='ARM TSI Percent Opaque') +display.day_night_background(dsname='ARM') +plt.legend() +plt.show() From a182317c628f873c12bafdf18885a445a93d48e6 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Fri, 11 Oct 2024 16:58:08 -0500 Subject: [PATCH 02/10] ENH: Updating tests --- tests/plotting/test_xsectiondisplay.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/tests/plotting/test_xsectiondisplay.py b/tests/plotting/test_xsectiondisplay.py index 8dede9047a..1c45db5e68 100644 --- a/tests/plotting/test_xsectiondisplay.py +++ b/tests/plotting/test_xsectiondisplay.py @@ -28,7 +28,7 @@ def test_xsection_errors(): display = XSectionDisplay(ds, figsize=(10, 8), subplot_shape=(1,)) with np.testing.assert_raises(RuntimeError): - display.plot_xsection(None, 'backscatter', x='time', cmap='HomeyerRainbow') + display.plot_xsection('backscatter', x='time', cmap='HomeyerRainbow') ds.close() matplotlib.pyplot.close(fig=display.fig) @@ -39,9 +39,7 @@ def test_xsection_plot(): visst_ds = act.io.arm.read_arm_netcdf(sample_files.EXAMPLE_CEIL1) xsection = XSectionDisplay(visst_ds, figsize=(10, 8)) - xsection.plot_xsection( - None, 'backscatter', x='time', y='range', cmap='coolwarm', vmin=0, vmax=320 - ) + xsection.plot_xsection('backscatter', x='time', y='range', cmap='coolwarm', vmin=0, vmax=320) visst_ds.close() try: @@ -63,7 +61,6 @@ def test_xsection_plot_map(): ): xsection = XSectionDisplay(radar_ds, figsize=(15, 8)) xsection.plot_xsection_map( - None, 'ir_temperature', vmin=220, vmax=300, From 7588ff78fa09500df0b461551dc24192a7e2c3a3 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Tue, 15 Oct 2024 12:22:10 -0500 Subject: [PATCH 03/10] ENH: Adding in averaging to the TSI data. --- examples/plotting/plot_satellite.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/plotting/plot_satellite.py b/examples/plotting/plot_satellite.py index 8303b594a9..c66c780212 100644 --- a/examples/plotting/plot_satellite.py +++ b/examples/plotting/plot_satellite.py @@ -29,6 +29,7 @@ files = DATASETS.fetch('enatsiskycoverC1.b1.20230307.082100.cdf') ds_tsi = act.io.read_arm_netcdf(files) ds_tsi = ds_tsi.where(ds_tsi.percent_opaque > 0) +ds_tsi = ds_tsi.resample(time='30min').mean() # Set coordinates to extra data for ENA ena_lat = 39.091600 From 9322e42721a3fc2d69836a5082d2c47630ffac97 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Tue, 29 Oct 2024 14:54:38 -0500 Subject: [PATCH 04/10] ADD: Adding new test --- tests/plotting/test_timeseriesdisplay.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/plotting/test_timeseriesdisplay.py b/tests/plotting/test_timeseriesdisplay.py index 508331f421..464dfb82f5 100644 --- a/tests/plotting/test_timeseriesdisplay.py +++ b/tests/plotting/test_timeseriesdisplay.py @@ -685,3 +685,16 @@ def test_xlim_correction_plot(): ds.close() return display.fig + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_plot_stripes(): + ds = act.io.read_arm_netcdf(sample_files.EXAMPLE_MET_WILDCARD) + ds = ds.resample(time='1H').mean() + print(ds) + reference_period = ['2019-01-01', '2019-10-02'] + + display = act.plotting.TimeSeriesDisplay(ds, figsize=(10, 2)) + display.plot_stripes('temp_mean', reference_period=reference_period) + + return display.fig From f175295c0b91ccca3b45f4293138ada4eb50ec92 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Tue, 29 Oct 2024 14:55:04 -0500 Subject: [PATCH 05/10] ADD: New plot --- tests/plotting/baseline/test_plot_stripes.png | Bin 0 -> 19671 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/plotting/baseline/test_plot_stripes.png diff --git a/tests/plotting/baseline/test_plot_stripes.png b/tests/plotting/baseline/test_plot_stripes.png new file mode 100644 index 0000000000000000000000000000000000000000..2943429c1f02678a837f5beae6623149d4be017f GIT binary patch literal 19671 zcmeHvWmHyO*X{-s5R{NoT2M(r6ai^z1wl$mS`=xd8>CS|Bn1Hhk&=||6bU6ny1Tpc z%>DR$@Av&VUq0u&B9GM%H#uh|qhW&}g!IO^(-Rc`A@7s;v@BV3Ej#MWx4CMNLQPzO_dZ&7xccoqf zB|fVp3xn!C($@VGb$g@!l>SRUj>|c6p2pw3J=^B=>~X?yKh1?JJ_s>!<-55#=t_@Ra8hdi!HPh4ZAuzJack5 zlaiC!x@nO&Z{LQ-#^RHZkVMDCOniS!70qYOP<617iMNV{T)cR3Y;G>#inyQN_QIct zfyj*F54WEiOf&>`^z=wsS>2@MF$-u4p<@&lrj-t+zGrV=)OBZ#^;ZNHuLXj9jEmF# zqkqK4u3bi&EFIz(5kbc9u<{1`0_E7RIDx)`7bGujY&y#wERWWk=sG$&9t~KuM&6kF zol3a1wRKak27wV73>KM}IB&Ba>@2??^*Gf|Ki2g)_SLC$z9%k@6&xJ=qoqYmUY=lL zVq$D!BFWyYHHx#z<7AW0*wi%HVRa<;+qZfcvDH`&MoDSu_*`f5O89!kY(7ECh@a7HaxpU{laIZYh(u6w&k}y3s zp%cR;rlW%kVc}GqwKA$;9epArBh!>qgZUg*KDHc4htLLsnFKspf9UZ}jNa^_VhYV(+uC-N`U8h1lWjS3rRYFhi(@kz}QoDt2V(l_p;-+94 zgvR}>=snUCngw7{$yPH`0-3-OCgT^2-GU1&IR6Mvy64)%6iEgM( zieq*#oe0M(Z*SkcygcL+v(MD5ly&1~uK6GjYMc5@$a{9se{5(tB&!3W6nuR`vg@tonDX$-dYDO+CnQ^%vIgT$Zu%eTL-T+Lb@7wc` z2ZIxv#9(V;;^Mxcp`mHo0x{?r(}E&Z?wL-eESTtKJHJ)+NkzbUD3p=!sFB( zL)3PmJH5>KL=GLbt$4X7Ap_pzRpYECHw)zkmbbY>$=^4%aLMeR^zZNrYZdB~J z^1Xya${sv>trmx_rlkc>$a%}NwDd0cH|vulr>{DdlyvTUiAl@j_5QG|KHyN26%vH4 z?Crn%;!}*Ru7-d^jvA1*x96#-sL*@*)Mt0Zg%B*Tuf&?MuD(9>&}`D6r7{r)M=Rc89Tnb{Pz087`?TzTA2qA>JmIo0zEyEp$aDw zKR>^q5m3(I=G>sq7vpk?(j52E)o<_4A|>VJY$qNA z1uwvxdVctDdBl0~hPt}?*wOyx-b}1zuFb4GwSfJ3DJiK3hwevvT$or`p%D?j@j@<1 zxw%BJ1kzwBxzF2a(R)7s_3@sJ3@(^}ugJ+B7M7z2pb`U>rh$dZqInaKty=4s6e0W7qZe zQq=nsHZ{Jc$g^8~jNXKTVN#eO;1k`2n?*H^EmR#&CPw;6crcuHkd|Wq9uaOxHp>y9O&CO zZ_b+yl@M!}+c!=$zFr4oPBU{j*m}_xCqO~NZ(9QOv5Cd3n}F?mKtJ$HqLP zqc2hO+u~L_Z!1x%J$m%6Ya0uJAm9}hMFFnPJB&#oNj8dulRJCjI8ET>-8&+=1Yzw? zzkrkJi~|}>h729eLX&H@c6N<%0uCxY!<8;->!+}#$T&DSFg;-jQkBYK%iWs}dv@2U zLbuxeLBccG2+GmX(YQ1;G>{HlkbH8IPK2zta#k}hoA%}EV|w1ab?c(8+jc}3|9+|< zJ~{iNpjHvGQ?T9eiL++n3x#wDS6!rAG_1d2p%}c6_?QQS&1S7A6 zAf}NJ(g~|H5Yq5}<`WSV6bu)33keCK;OFPB8y?nvp$p+*wdzp*#W0#AAh;A16cGv& ztYX6MH*VapUK_h7c{Kk!m4TjKl&k0c?Ea5X#kE zX8lTY^Yeia!6gM@5c8CgvNAG65SNqm>-|_4r);=f_pDL|<2bd;u0ozGku3F=)pp6r z$qDt32+_#Mh(Uk8(FE9*QT*m~99yGJq3ygAB-_rnl< zLwP4d%q?nWeX=P@)I!s!E`CBnLJZN5!fu+&MGs?n$dCq+(<9}Dlh(sC$4zt| zDi!t-S5ILe(RWNE8BY42C#3(38}g@fXGQ|9$wr3rXKx7KYj15OL3|Pu69={-MI}Rg ze0?R;-EeCaIs#8_sDC0+OewG!)|^{hWOMVEj{_$Gv0-#n*J${6s)BE9tXBHjGMKI| z5%%e+Q{0$aMnMeMPo-1Kv=||0WJYjtaBw14A3)-YhZxLtBnGKVN%?hjbzL+cEK+rm zL!&3s@I1loOP7b!S!WH8mh!CPTW6giD~;G0E;84$K0erfSZZ?1*xK4^_@_hS zFygkW9bsqM^Oa>+kTafSzceHJ5%I)+?70PtNwwMJ5aEl6gn>bl618UsxWnC*Vdehs zO-<6sqjLE;dLNIY2Lsv&rjySuv^%y*>(t&X>-7KmDc1*E%0s(y-Muk$6Efb0^@ z4v|N$S!0qqymAj8`g(hNcTq`2tj~2aqG43u*f?#caPKu|6%o?VovuvAweRZa;*wdw z$<9tXFffn=Ystw?TdE*(tb!0-xpIYd@rkx}D|@zp{Squ4WJVEXryF=%P;n(@W!-=m zRH?GKWjP|)=(JsZdaOVkRo$tmX?-K16{?L$SY^B|h<3MLFG8+`AqvGq36s181gUihx#^XK`P>jxdB^wGCp$VjYoR>) zlMPeTFkJ2s`r$+RR%R>YcqW!ES-RE26O)s&LLTsxpn?MKzP>)@tVKwJT9$J?&D;rS z1q0v}_Mnv6Swq6zIHC%(_ex#C6(jjtIe+3>!Qjt+a^|2DawdfZZh#?%#>WT4ETwEW zH#Zx1eDYF)IC}v+vEqFFuV3<9*mjGuCOuy-L(<`0eR@KLLqJ^vh2UElK0QA_WoKun zdXd@Xo#jEuq?dwBMyuRX23pfw6CxwYAzelNYy)JcBO@oiI=qm{tPeYkCcgOnWp1=u z1cGY^?ju>*N zU?Tj#eY@>`FwfK$a|a*NsaHisJZQ@1dbnI{GAhR^KjQGB>6y8CR(Ezu{)TG*lja+c zm%+P;!!K(n=p&k<$ySRCOeaSwU5tWrl7vr>Y`9-^oq@nv1L>@yx_SUS*#O)b)V$Qc z5=02`@n7@8-aL5Mq`T7>nD*?E+hIaL8nrhL;diK;-+lNHNF!jcx3@M9?e1it{mdM8xZ+3@9zp_0@M?MkjqYVrO1+v zWM{m5L<^N6#GddEbXU$SEidaX@_ptPdAV@{X`cFwAe0>s@_v54ANDCT^EEj3Nm%N- z3a8DR<^$YvA8tv`<8+P`z;fxYPc$N)$A`OW&W9AoD)~uw-9vjXK-6^w)C7Go7@Rl6 zAQH%A>LPDwp()VIS+V;Y!qbK3}f z#5aZvR-EadP%1O!fBVL^{j7e`YU&3uln!#&PReFelLUo?Fg_{GZfr01c7DE35}~{9 zzO(dUU~%u^`%wnY>fuX{$Gsz9rK$#ypQ*O_3N~?Wlqr7n^*xKGoNuXk(M$mli+V`K zTZg_t76QSj3C8(nB$C79wy3-G@5NiDJ?uxj)u+DYN;sq50;cc-HK*10May2N^q*PSGC+O&2LVC}({TokoCu9s< zC)=?h3m|WRrJ$P!aJ4N~0!d$ZF>rL&kT#$3Goi8U5rCP7B~=H0~>!_GGhGXfeK zB+zVhxePKuVlFN&z{$inXMRTXXqUcB>0bShRZlgCt;1B;!^0qJl(~PucDgl+n1sX|_6wW4 zXcpLeZeCuk-vz2q$;kkAQHlbx3b}rEmw6AopNWU)0}JpSH+gulk*80ej(uwvR{q-h z;|GF3UHHVhe8uM+Nk^I@Ba%-ywK!bP%WuDg%gf7)#>?BcAxeL&0f<4**f;>-0F7z? z(5Zy*NOn@smRUUjb%GSgb9Ac-M0d{H?-NGd0UQ?UMrOkqqw^OobXK{!pd}iN=@mea zkYaWuOA{8rn?nIu3qf^#e>R~WFg?hmt!!;2baZGhl8~Uq2%6a)hcBq&5fFS&O>htT z@k0hJJJzXf0gw9%ILxxs@t5O+#hkgJ(%XK7bitCZ$rAyYfp`4AJXnlECQm9Hi2)^1 zb6vP=F~kS80A7|>u@nFV=p=#dO~BjSgM>8?JYz3U4?vw20GaH~C5Mdt1wlj=6bPZ> zg`$C?tE&qIr&ZoWI%HRT1(@=aT8_@#pG=A;Po7w<4BcHCDkXuspVji0%8Xs|r%wn5 z28NZDl_;d-0OW=E6Ve&1j#QSo?%zVu1OQmV1**UlC&6meD;&81lZ8|v2h0SbM4sci zL7~Xj9_v?)Z`iO0{2&kFw_EUmH%u?Ewz8^$jEAfhaLmSy=~fiD-rL{5!p?pkp5Iqy z$Cj;A8C2zVXi(#c!FtG2MYvBDG}#oQHl)gB(14pl#~87@E6)%L3 zs27?<`J|KgRtRT8HWeCq+2>N$D$?0=7m7 z2@x|>v6fXVyFZB;GM?$IhBS~UKHRJ)6bUFEb=TtK=MRD;4ef7?udHy)w0%h4v3`I8 z>@dl70FhjxqzpnqPe}1u;%~U@ENP_dAU-)cIWmC#j*pKg!D}W#9g?J)`S@S#caiJ) zw2CU%L;EoN(n%CU2W|>i1R@ZJrt+!4ncUo2p}Bk!NCFAIH@?^8JqyrEvQK* z-G$xYMX4&UM@lK&1;x4>!2l&8)c68Qj8>&krzQcIY|!&n!$eoNgY#Q( z!HH)-V*;V@eLSrg{c?I+dBeHC%B>uQtkFgQiivvuF-k{PS{0Dy z+({(%Q<_OaBIg5OQ1bDDXhS7M3+E0u_edXOa{wzFlshLUCxGxUZwBLLR6Wkod5gld zF#AdepWg=r%sosya9ImT`|i9phVk;{OO=UCwVdyezd=rT)HS*YjYS5A3>Gv5Ih%t` zb!4crp-5z7<5@^LP)IG9`YtBy6sT@;GW=yE;B)mdRIPJ`c?AUx5JX=C7%@3k>y4I5 z(9MIG;RlIx)%)4ui}zKbDE6(&&^I%C1COlQStlYU4&7Q8n1j3I)P}Pt--M#>LCb&(>ctFNZ1h{i_Zv8?k6x{-0d%(jPT;GxKHXEoB3CK)Z(KO!g(m*j)}!6 zfQY)U+<3yjb7N@Wjh8zO6;k7U*j%!+%fNSs?yL|N73owe{Yf2RuBG5bXY6eZgm6={lkF8x&+fCx zloXr~{&|v2`nq~Axurx~(=A}YuPC3#skwxU)0 zR_HzzCH~x8{>7m~9enFs6>26zxH9~GIt;&~VUy%icSPgH>hO0FrYgzXISTrZdI|PRi>=oC1i5#p z@xA%&gLtA8eo4!YuHWR2G>?k1!Jj4gCF$RBKh4x4#N|UEUzH)htFvTH-_ZL0euWd> z(V?7Ipf6+o>UH{s5Zp6c3D_j~{dW`M;yvA{m>6><%vQfQ&mH{Yj@IO3U=zY9?^O_6 ztC_!9_-EJiv9HE=J@cr{j+3upy+3zz1KbmDkx4y8<4LPvB zStR>;Kd>~cJ+`cLsmp^>?zMTb<>b>dzCsDDLR8nU$)#$PxA=4yVIfz^su8GC?4HskX&)wZ}ihTqBqb+o=Jmwq zitLOxXtf@vkID~6i($+Ar$;1DtYX@2YUMx7(a+~<*{z&Uz3=F-MXKHt)wuHc>ewtD z!&+_Mx!k-$>5b*9>>4-mC!1UmEv%V3-)DwaVp|u8TV_ckIL_Oe278UsNP{SoPKk zpIm>dN%vEgIo9VP`8q=U=LUPH5B&SH65VG-NUF{%M&oo^S7R=GFD#hkBy zIVmF{u0N<7$)H8(xt2M;oNVD#QQx_=!}C~1?tZtfkfU&r>RGdnZh; zJoe|8g5qpO+0{P>G)J8tu1h~KuhJ>GMDqOzXNPAyH2mW$jo0oC9wH~W7lIPX@hJ=^ za8_%psv^wfSYAZuJteE67CN)DJ#DDz9Y}ZHyYj7>Twrs0blwzGW6AWO+jhNX#PanJ zSDc;g8&A!Zesh`kJ4;lMY`DtzNKdpr;_Qp1Y94qxGm%Ln%{7urw)FObk|9*c1 zu93{MOK;9wv+)wz1jjcqWl`4sdZcU5?nJm|zNA`VY3hC8T=;HdjW}{x{&6gmT^`rl zI6c!xRK=V4xpNy2U5DDv7GdWLa9hq&=#4in?v*$N_W3UERZ9`%rr!{7s49=~B-EzX z%gn(+WZc(?Ugt2SebhCK;je#uy6;YYY$M)T&0xR(B`GfHZOrm_9tW?iL}80rp_*6t zz3Zd)5A~m`a5L?0tZz4abjL@9xxlQq;Ras58(p+bbUg$n!Y)tn#;E80RYLX_qfS%J zYEiy!Hs|BASl{&O)yCPAtcdJX4aMgw-?zMMu==_=V!~eetG{GC$H<`HJ@{pP61zUjKo+8PZm;i#AK(DgIO*bMEP$PujZ8QSxRl?Cma8JdxMn*^&hDW&W) zecZSY2YhSOUl1WRiT8;poz9CghG1A3GraCEJzLXQN%v&vu=dvt1ks!9x-9?mu;zQp zZKP(PlWynjsGr+pvHtQLA(ETBc!ny)=&b%8VWLuHIeO?L(s`|>Nj z7Zhid*%UO|x=6@*ylxgg(#bZK%QG84hhI*w-gmaPyGuD;#7XyA4~_l{y7RJ!deLNW zo8kJe)UWha%|2dh(r*7kV)ju?f08uFnP=od-uu=@!Vl=**>Y;qVIP_lO~FU0dJ-xV!!A@tG8G{}?38XuW=`l`4{ zRIqb7Kq~9?jY9$nQ&^gAeu?F(xR~bbxa-SxW8GiBnl|dr^HFl2_cYCPw=(}jc!3g| z*_>%q%KhtuUj1aGq;LVJ`0;VaCq>2b>n%1(7IF*OE>YyQ#v|)AB?a$ploy?slFjil zp1(t6acxHa#mIv;qxz)Y=!u-$Z#Dftsjiork-d5xhn=LF>Se6LfD^Q@_FF1GY{KaU z7>RjM<@e^O*G(~wJ~j20^nyp&NrC<^<*J6d8PAv><0NIBigMq8q~I$xzNFTvRlzo) zv~tL!lxTxebLN44|KvCizL=tj#+c#@lI6d~jUN{bebr_bEcdDZCd(I7%x3bNgESsv zeC<3lFHdK?mk9|`G12A{+(kNN8cKeEt)<5Dl9_Ie*2Ls4BX{ToRKPR1UV5N znfon7>Q?5T^etp*n;|V32m8$775IUCG98@3-Q1Oku4XDLGUs?*Xz}J|0uS!y3-9bK zs#XKqpTskxwAw;X&5zvNIx?ksw-!1toK)raUgjJ{e75FrICcB)jP=cJjpo*V2^dxD~Ik>)4FxHJw{#epQO&R(NCART^t(Q`v0}cHJhW!K&}<`o=uE zEZ%vAo|*=;uCt;NrkH*7%bq=r>noyebL3Uo1<%xwB}XiNr$mKfJ=^)KhXR?Vv0rdE zE&Xqv+RrU?s`Px%!Z;|jm8_1;#)OIG|02;i}CYE zYUIB&y*aSI8>?wOe&g74N93`^^GGlqedfX62mjGvB}AS*6B`46q|Xsq>Hh8 z)>xdCe>M$J4=-nR4v$a@_vAY}rG?NEY>%7Xa0yEZ^W7-@-CReS$PBbrncLXkD7S$D@u9Aq8-?LoQlQNcBR^PhlkD{$CT68 z{E=UduK4+o*u>oHktx1oZ0k(Jw_p>30TZe^-6mY~LLuXi#3>WQ0&1+1{2WXzhAwN;M@MK znHbLLdCC9%=2Z+R(*Az+s^|Z|>Ho!Z>S?b^=*c<@h!((8Qupp*K(n+nS6}?n<;$wu zxqnSy?pmH`4z~P9KCcwnCTkAXm;nOYu{B`Z)&&hx3%&xByrf<^8oEm4|tv zQMuoqw;KS?%Lm;IVAi|=6cD950#nB`KAtu|Ki_I|=HW*^YfON28=!5q)<~f{7N#KV zm7Q&&N|^BX(X3u)%X|tgM`-~#+XFzl69CIwgEO0+fQ>qCd!rmt z28N{GAb8}tD=6s?Fx$LtwB1kx$jZoA+CL8nn6uOK;?RdaK46+ykILA4UTFGZ|S~w&GzR z4m|?+j`!xj``KAr{gZ4oG?;iSN1j-e|BU6WA09s5Cj9qm9V(b99_Ztm_2ptfr+}^S zXS>Mh8|XsoL7(-1jPQlpHj-uwsCgaC{!O^x2t~U4zcE2W3Il;o(B5KBb*{^5JlG0EnQUq=mq< zWetGS;mby{-;+F!KD5Rv>~p% z&fbxEZl7s(!Pvy44oG5v5jlDPT@qN6K)3|RT8~zC+n}?9;{BTx!+$+p#uf6pKMg>K zz6qRKzz35(PET@yE@WcD>VA6c01Zq%;G_1HJ8%L@g$-9%Z2l5Z6?*#5U>Gvaz||5; zN*7~z%!2^gXLWl|OHAwqm=Yr!8$O`G_W2c!vU*7o7n`h|Gdj>Ci>5fTwi}?yA zCFNeb$B8#URhS6MEP&NeX(akPy18q)qB+-r+xZCkTK2;Ln*lfX0&?^AZEvE9hXil` zuBi?NjKT+2{MoOM>C6q&u*+^5em@UHMl8f)xa=$&TPEwJD}QaP^1M4vWy~HdB@j5z z0MP&_wU?oxQ(`kqKtRps$)pgk2Q7NU6Cg}mXcq)Y%l_Rv5sjFYF>5cEJ4&{a{Lhv} zEy{X-!w3je>DtdtOoH0l+IR*34vE1z0>(rFko1qAKUK>)?%ur&P_w1u-(g5>v4WGl z4kQhgDi@Ao=WP?{gH-O{`}_KxZ1uBnE!sXQF6ISi>kr*ObQ_@LM=*2fe)ym)obCUt z_8Q5S3D^J@@PdD|w$=k^k9Xw^tZE?4i9Ue;(0SKTB@F#{*)c`#$xP=D>;@n zIMv5HT4cfD|BjIy?)a^^SncQ(?eAr(akeFh=n?}@@rtXP8`@AnTQuxl|87Pu;b&Pf zphai^Dr$j&PDt!6SRGb3Ioc4A}K%bGXHnzD#z$f0TMSCSR%w9 zaTbU?_Gkxo4kn)L?VEqca~3s!^%AQq(6M-h_R!?x_+h^{fETqoJ#mBZVF1(_RUwIg zE%F@dSWHb#QK}Qb;8M1>+(1X@1jrkD$cDgB2L3h#)Bh|93{V1EcEC2QgWhe1GxOi! zZ3O{>eGyjgw%JHUY27k#6u6$XU22OLI=mAc^zYbkz%%!NrJtu2yaY_Gn%M-8U;wfw z!JP!@nsG8S*Ptw4gGM~qm^3>(J9}XN&|MB;rM>&t;8e*it)p^*jbS+TBUC+C599Gy zb+5ZNGiWOo{0y|P@9FDbK}2!A;H{xydtnZ4PZkpM@AJ>#+Xj9sD*6uZ`p*uBmM@?C z5ucKh61X1(nn0+)1sHU>g@suQi5LHQW~Z#4m(D|I64gNO4M>&$y_YYb4T8U)|DR(1 zKb>KBWG`TJ#HxY~ApxheX$+BeOwbQmUg8*fSoLbD{2by*MPYx1aT^bzXA&oWU0s0l z;e}dg><3t#+{>qcH$Pio5a9I)IU56YDWaT*%2J$|KPSiq*dgG-92BI{XY{j|A zkI;*#k;_)wU2?k13oVfqCM@KZatOnz4kq4)@>9cy0>L=DM#q9}3O6WE&k6pJZK%g- zL0crs?5PQkRhCX)(*S?^klwU57H`HeW!NGk;_A9KQ!8n4axi|S;EHMreWWzrGerYH? zA*z>iNlYqZ=}AOFyY+L%(9n&ekJaRC<=63n#O+X-YPYi**CH|K*@SeKt1NZG zwFQuOv5Gx~7gKDru4)haFWFA>Y%Bv{ZFn|F+g~lnJW92Dg>L$;BJXv(-sA1f@^_&x zr}r!327Sv}ByE%oR$5|;GoQb1T#39Dl@}_b5ZRt`;Ff8c{1HoIX=60B#O5$);*NaP zc|nii!&%#_D?{%aZ^=SqC`w+MAxnbF>FW}YK&bs<{s7F)7xh?fjU_Qf_GIJziS#&` zhwlM9=LjceGZnBiEOXU#iU*9nnP>LT(&1L8%urC2xdW4A(Qv5uA(00oL36RutX_8_ zV?a5JsmEf`QlWl(-tUqaEBrtf2LvfMCSKXWU4Br;1uru}`ReCxoq zq0quWen;G1%OUoG;-$jzxmve?YU>LulE(>r6*z&sMzgLa`AZ)~?4A<7pD@bfc3N8F zxZ0>IshPjTcOqOSJKbjxBFS*~9=R+haTsTU0C8Zrghp(6dj6!ai{QoPKsD0Vg@Lr# zQ5z^j;+NEYiom_Iq}u*W=Mah;w;V`sNzrJ8emA&a}; zb+kFWur$)g)ZrTja=(8U1My#dIyQ6cdR64z@!&+a_n{1{D|2jXlW%@s{h^HW5$+3x zr>#sY(+um0Umu!`Gxk?aT8s}2Q>VDC3U3l-k5qa4*lL`*0Kxmu)&9MUzgJoeqz!|9 z#T6t?@_CuA=R;dLGEQ6(hmb;m^L>QDwW!y{{SJH4T_36W9MAMxGAER~gv8esI`)@d zqw^O0qAvV$!znh2tGP+L$KEA4q4tkNzthL)m~3}>j^TzsNk>&fDh|=4U}y6_fbxMuL>)UtWXz`8H#Z zOt39RT5%A`CA|ghp^~UnQh|@zt}E2??P9K->hI}(956e%(>zcfF0Sy1=&N*%FiPRN zmO3eAL2E&=eL3Qt-(vsnQyZo3$HPvgvlADdJp3tY;KyDy>thoCH0MQeLiLxp%2%1q z`F7_VzCq_w+-dY$lkHVT{Ple>Z~E9Ul3i6BqRp#>#qmO71K)tz{JVX_T2|b?a+o+I zqA1LgaTVV)NyO=`sH%cE?v)H@qOF;)TEzrV2W2ei|I1(+pU9 zUceMnYPxAV_Xq?~C=&^EQ#C+KWOF^A66(AHd5J#|acY1}4x9y$KS}^0gX{S(YC(4A z&r%xE<-N&Cy~98jN7E$`;ev{f5VSLnr^maY(b4J_ba(DJ{JCP;NJdIp17vNq5(1G* z%X0FUnG0#OrJI|Z-avF|Nf;X&dwHn0q7%~c!~x2MXhEl2kmEkee0&4g$+bX>1VW(E zj-s*+({{zRt5?%1M_8SQ|8!SibhXnUsCEa$CU;`)m|}w%4&{?wy?O>@No>U-F&_n* zFRHMhyc0x1POgn=c^8yzs7#}-r{^)c$_Bt?MU{eLT3T8X($Z{&G(2X}pr#}PN)iQ3 z6mi7;O+X2J;^N{!@SFs*>W^ytKwSN?qr01QwQ@K9!69&+fzH4l$)U{(B^Eshf3F|` z0E(~_#)9f17vwjnIs&9xB!~p$*DWVNy6_1L3xhn13>_PQ4Hw1z&J!aIKfw9C06f0IJ zts8bnJvhH+W%+`Z^~aAN%-<4XV#*iZ6;Gg20@NrlL|u=!d#87ll{c?TZ2yP@7L(P> zmlshU+kp*ma$g?mujo{^OsxP76{;2S?i^iRcFh7An4s6F8W`x8uK@;%>{odBGp>7c z2dIz{rpurSFgV&<$Mgj4*^3Xy-Q5-Et{e}YyMhuz*g&Pw2!_%KcWhj(KIH)U9pk@( z-`1tn>$CD1P=zU}YWo0ZzY?Jb%F&kqA@?J7AZkk#nD)x2F9Y|P?DQ10i;puz+gn;L z0)Lg&dfy%RDWrKo$Hqk_Ku4r%LUrj9GpHRnQGv3(92kaOI>k7U+|D2<(>28r=-Mbb zslUIUOgQrW=lF7v`-DbCrRG<Jk(92n_2mjDYqhgvS+BTr$=uwRH8WVZep^M zAwJT*4r>o3CIc^Tahv^}JCy{Svq1uUd}DDABZ{S2s_uufAU~H1p}he4I&6J9l&iv!Jep3PCZj*tC}b0eQK33fE;yN|XoozOvo4qQTzBm?C*v#v>rHhuaRaQd;au~`ROXHTFUH`=KmQeC!7 zl?P%u>0tQu=zD?mnJiFBLite-iYD-*A$WA{)z>UNZEbdwqRTT{hk;3-3uO+dpHY7i zAsp`^yYe}b*AC{)3;r(&FdCrAS$sG?(%Qh5;dB5A0qiy(sG7TWU|;!u`Engd&n-){ zqu~PSdR`~*K`+?>@Fu%@E`sVrp}=_;6B7g*7p2yv9D^PZ&-wi|~|T8IUDY4q4dQCBSW>Uk08n zTV;EfO4JYypA~|16ukHyL>)P>{GpsTsvHoOPe9HDt??xRdvq-SD+XKro|gke;z3B)PK+T^Dzp1ld>U~W|46~QSr zf__5ZzyK3aDxiqzfsf;K{dnhwgHmfOKX9l?e26Y_8Ax0UmeY9HN=08eG(R^t=1nME z1FmrmHpe&!AY^1@YCv@dN_3x4M(L^!>z;<8`;z^1i;qfC4kgsotdRg|1wlxu1=d~l z$pIS>G?Tl9`S>ouOewKHH!>RcI6b5TKtRvPNYz9H_<1)$J`mVy2!G;qVNMzN_MFiD zWdAuTJOCg|nKG=*p<4Hu+`|~|TDU&|e?Bn#%oD(+a~gMx{K7xA0wh2b^#!qL%8q1- zU|U>UY>T$0CM7sAvQ}W)N`uJ*^N`KKGwnlZn8658LSL8K%%X=rz-F$2B3E~B&okx& zuV_r^CO}w#;G$Gx0QZ0geG!(t5nh&vj7)XDyg6FnLvNhGeyX$0Oxrz(XLHyx=8 zbf`QTfo+UM4b}jtMTGG8m>7ZdF{U|A7ND2|GuHq0YkY8bK-jzh%FL7lAr~GL z&_Y$zTOG(vhsMg4q2||EV2l9%CV>{+aWShHW#Ox9X{jB^N1{_Oo%I>ulPdtK0kp~K zaJm3d0i&yJ7~XR>JCr&#Of`KcFT6=)kwkHw{I7&5*Y!(-*dpBL+3HXV-RkGJUB|zXK1*IDV-D_ZW9H7IILcr^`bU~_*e2@BQbI?`onedL8YF8PVsRNR4Fiy~DHuI-BPp(`N*Lk^ zD9S$I(P=&W^^XBUK8CJ7XeZR>9D_2Oov=Pfm7}Wejj3jkWw8R$`wU2|0ma6G?uBPw z9v3j9!_e#7<5bqjccb3H!3nU5(DPz|HUFpqM0KC}h;<`Ci9cEAi_dE91o81%PC+$5 z!e`E$sTmoe1!P-sU;o^GbWe_M2skumx85nECQ^`?0|gn?paT)qD>jyj$6^Q<0a>C` zqFdzE>7i@auh&7cD25A}LVZt9kO9S1meCJiKkkSRY-JTZB>@;SQrH3cVx zLkghVvYInJjPZfv8?Ny1kgu$)811N`Q+3;6I(?3b6~eJ{=TOA9-Hj_}Y;a!8WR+}f z*`g|U^y~sq#HJklGbJC&VVmE>^#Y&ITKF4)Fi4t!<6n!nt@hi7>Dtkm<7>vo#-g~1 z8!MebWNfTQ>9Ug0eZhxwR1^=k>){>8qIBKw0g)Ioq+_+7UO*vHkp<+1AaA39*qIWH zZe~!irz*L^!2{^~3=~_Q))_)h?g_#Jv>OKLHwrIWIXHM2Q6#F^L_wDW)&`DiAO%69 zXJw^Z%7MYkXS-o-n8RYr(d>@haua8CI?j8<5K3P?L3qi6uh3EI7ZJZc>DDr}G` z$D`Q`oWx+!yAzI0PI?2m8}+^*ziWcV5u~T5ena=={|LHQ9vnavIEz@r z*&1-34j21uOP?G@sW#Hpehd!rb-{9bDXihR0X@P(Op@(1)zC-cPSwd-PX)h1u zWo1d=LF6FoKr0JJXXhYrTNwpmyD2LQ>s;x|Jw-@EhH7T`z-YCHF3MI{*U&(Zl7VT* z&vQR^jyP@!Q68q#)zwWY5^~ug?dj=3x$W?+Kxodw$x}WrXHx~0+O;)o&WJ+U9Nv#( zqh?nMYfJ(a_r+hoeu=@gAbiGSVq)Sl>1GZM4_Eo8a8zoK<%cFD94-XCtX6X^s8xKn zY~VlvCg`yEzJ2>PO*J#J1A-*vy%eCJA_7;-UQ9TXy^vFmYT@Cw=*ccnOG`pyr4fw3 z6)ilo%chngu))bFf=(|}f^)viJPcR}Er%mQ(8fY1gb4I1c6N3jI8~&t)aH4qTldt^ z{RlrIF;CVdFz@dWjY7l1lD8(WUv+EK-iF<%(58j{(1pJ>RKO3KN|TDl%2y{Rhn8>- r${2Kc#6SV!i1C>5fA@WkQ$lOLEbhfE`E@vC47o2UFOe;#_v(KE9wO<^ literal 0 HcmV?d00001 From 7bf846b19c307549ebd5a01ac76b74a65db6f342 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Tue, 29 Oct 2024 15:00:11 -0500 Subject: [PATCH 06/10] ENH: Adding climate stripes plot --- act/plotting/timeseriesdisplay.py | 120 +++++++++++++++++++++++++++++- 1 file changed, 118 insertions(+), 2 deletions(-) diff --git a/act/plotting/timeseriesdisplay.py b/act/plotting/timeseriesdisplay.py index 45f321aa08..4b6c3ac0c5 100644 --- a/act/plotting/timeseriesdisplay.py +++ b/act/plotting/timeseriesdisplay.py @@ -9,11 +9,14 @@ from copy import deepcopy from re import search, search as re_search +import numpy as np +import pandas as pd +from scipy import stats import matplotlib as mpl import matplotlib.dates as mdates import matplotlib.pyplot as plt -import numpy as np -import pandas as pd +from matplotlib.patches import Rectangle +from matplotlib.collections import PatchCollection from matplotlib import colors as mplcolors from mpl_toolkits.axes_grid1 import make_axes_locatable from scipy.interpolate import NearestNDInterpolator @@ -1847,3 +1850,116 @@ def fill_between( ax.set_title(set_title) self.axes[subplot_index] = ax return self.axes[subplot_index] + + def plot_stripes( + self, + field, + dsname=None, + subplot_index=(0,), + set_title=None, + reference_period=None, + cmap='coolwarm', + **kwargs, + ): + """ + Makes a fill_between plot, based on matplotlib + + Parameters + ---------- + field : str + The name of the field to plot. + dsname : None or str + If there is more than one datastream in the display object the + name of the datastream needs to be specified. If set to None and + there is only one datastream ACT will use the sole datastream + in the object. + subplot_index : 1 or 2D tuple, list, or array + The index of the subplot to set the x range of. + set_title : str + The title for the plot. + reference_period : list + List of a start and end date for a reference period ['2020-01-01', '2020-04-01'] + If this is set, the plot will subtract the mean of the reference period from the + field to create an anomaly calculation. + cmap : string + Colormap to use for plotting. Defaults to coolwarm + **kwargs : keyword arguments + The keyword arguments for :func:`plt.plot` (1D timeseries) or + :func:`plt.pcolormesh` (2D timeseries). + + Returns + ------- + ax : matplotlib axis handle + The matplotlib axis handle of the plot. + + """ + if dsname is None and len(self._ds.keys()) > 1: + raise ValueError( + 'You must choose a datastream when there are 2 ' + 'or more datasets in the TimeSeriesDisplay ' + 'object.' + ) + elif dsname is None: + dsname = list(self._ds.keys())[0] + + # Get data and dimensions + data = self._ds[dsname][field] + dim = list(self._ds[dsname][field].dims) + xdata = self._ds[dsname][dim[0]] + + if 'units' in data.attrs: + ytitle = ''.join(['(', data.attrs['units'], ')']) + else: + ytitle = field + + delta = 1 + start = int(mdates.date2num(xdata.values[0])) + end = int(mdates.date2num(xdata.values[-1])) + delta = stats.mode(xdata.diff('time').values)[0] / np.timedelta64(1, 'D') + if reference_period is not None: + reference = data.sel(time=slice(reference_period[0], reference_period[1])).mean('time') + data.values = data.values - reference.values + + # Get the current plotting axis, add day/night background and plot data + if self.fig is None: + self.fig = plt.figure() + + if self.axes is None: + self.axes = np.array([plt.axes()]) + self.fig.add_axes(self.axes[0]) + + # Set ax to appropriate axis + ax = self.axes[subplot_index] + + col = PatchCollection([Rectangle((y, 0), 1, 1) for y in np.arange(start, end + 1, delta)]) + + col.set_array(data) + col.set_cmap(cmap) + col.set_clim(np.nanmin(data), np.nanmax(data)) + ax.add_collection(col) + + locator = mdates.AutoDateLocator(minticks=3) + formatter = mdates.AutoDateFormatter(locator) + ax.xaxis.set_major_locator(locator) + ax.xaxis.set_major_formatter(formatter) + + ax.set_ylim(0, 1) + ax.set_yticks([]) + ax.set_xlim(start, end + 1) + + # Set YTitle + ax.set_ylabel(ytitle) + + # Set Title + if set_title is None: + set_title = ' '.join( + [ + dsname, + field, + 'Stripes on', + dt_utils.numpy_to_arm_date(self._ds[dsname].time.values[0]), + ] + ) + ax.set_title(set_title) + self.axes[subplot_index] = ax + return self.axes[subplot_index] From 9b777e337f055d224fe193f44fbd4829d43b845c Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Tue, 29 Oct 2024 15:34:00 -0500 Subject: [PATCH 07/10] ENH: Updates for stripe plot --- act/plotting/timeseriesdisplay.py | 40 +++++++++++++++++++++---------- 1 file changed, 27 insertions(+), 13 deletions(-) diff --git a/act/plotting/timeseriesdisplay.py b/act/plotting/timeseriesdisplay.py index 4b6c3ac0c5..c08e712cde 100644 --- a/act/plotting/timeseriesdisplay.py +++ b/act/plotting/timeseriesdisplay.py @@ -1858,7 +1858,9 @@ def plot_stripes( subplot_index=(0,), set_title=None, reference_period=None, - cmap='coolwarm', + cmap='bwr', + cbar_label=None, + colorbar=True, **kwargs, ): """ @@ -1882,7 +1884,11 @@ def plot_stripes( If this is set, the plot will subtract the mean of the reference period from the field to create an anomaly calculation. cmap : string - Colormap to use for plotting. Defaults to coolwarm + Colormap to use for plotting. Defaults to bwr + cbar_label : str + Option to overwrite default colorbar label. + colorbar : boolean + Option to not plot the colorbar. Default is to plot it **kwargs : keyword arguments The keyword arguments for :func:`plt.plot` (1D timeseries) or :func:`plt.pcolormesh` (2D timeseries). @@ -1907,15 +1913,11 @@ def plot_stripes( dim = list(self._ds[dsname][field].dims) xdata = self._ds[dsname][dim[0]] - if 'units' in data.attrs: - ytitle = ''.join(['(', data.attrs['units'], ')']) - else: - ytitle = field - - delta = 1 start = int(mdates.date2num(xdata.values[0])) end = int(mdates.date2num(xdata.values[-1])) delta = stats.mode(xdata.diff('time').values)[0] / np.timedelta64(1, 'D') + + # Calculate mean for reference period and subtract from the data if reference_period is not None: reference = data.sel(time=slice(reference_period[0], reference_period[1])).mean('time') data.values = data.values - reference.values @@ -1931,8 +1933,10 @@ def plot_stripes( # Set ax to appropriate axis ax = self.axes[subplot_index] - col = PatchCollection([Rectangle((y, 0), 1, 1) for y in np.arange(start, end + 1, delta)]) - + # Plot up data using rectangles + col = PatchCollection( + [Rectangle((y, 0), delta, 1) for y in np.arange(start, end + 1, delta)] + ) col.set_array(data) col.set_cmap(cmap) col.set_clim(np.nanmin(data), np.nanmax(data)) @@ -1947,9 +1951,6 @@ def plot_stripes( ax.set_yticks([]) ax.set_xlim(start, end + 1) - # Set YTitle - ax.set_ylabel(ytitle) - # Set Title if set_title is None: set_title = ' '.join( @@ -1961,5 +1962,18 @@ def plot_stripes( ] ) ax.set_title(set_title) + + # Set Colorbar + if colorbar: + if 'units' in data.attrs: + ytitle = ''.join(['(', data.attrs['units'], ')']) + else: + ytitle = field + if cbar_label is None: + cbar_title = ytitle + else: + cbar_title = ''.join(['(', cbar_label, ')']) + self.add_colorbar(col, title=cbar_title, subplot_index=subplot_index) + self.axes[subplot_index] = ax return self.axes[subplot_index] From 48f80295f1a9c3b2a7172f845ed3e2fb656fe4f4 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Tue, 29 Oct 2024 15:34:20 -0500 Subject: [PATCH 08/10] ENH: Adding example --- examples/plotting/plot_stripes.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 examples/plotting/plot_stripes.py diff --git a/examples/plotting/plot_stripes.py b/examples/plotting/plot_stripes.py new file mode 100644 index 0000000000..c6def32e80 --- /dev/null +++ b/examples/plotting/plot_stripes.py @@ -0,0 +1,28 @@ +""" +Example plot using stripes +-------------------------- + +Plot up climate stripes plots from already +existing climatologies from ARM data. +Author: Adam Theisen + +""" + +import act +import matplotlib.pyplot as plt + +# SGP E13 MET data has already been processed to yearly averages, +# removing data flagged by embedded qc and DQRs +url = 'https://raw.githubusercontent.com/AdamTheisen/ARM-Climatologies/refs/heads/main/results/sgpmetE13.b1_temp_mean_Y.csv' +col_names = ['time', 'temperature', 'count'] +ds = act.io.read_csv(url, column_names=col_names, index_col=0, parse_dates=True) + +# Drop years with less than 500000 samples +ds = ds.where(ds['count'] > 500000) + +# Create plot display +display = act.plotting.TimeSeriesDisplay(ds, figsize=(10, 2)) +reference = ['2003-01-01', '2013-01-01'] +display.plot_stripes('temperature', reference_period=reference) + +plt.show() From 696d20a24e81525c3b4fc7f47ce2d904f9263b49 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Tue, 29 Oct 2024 15:35:32 -0500 Subject: [PATCH 09/10] ENH: Updating documentation --- act/plotting/timeseriesdisplay.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/act/plotting/timeseriesdisplay.py b/act/plotting/timeseriesdisplay.py index c08e712cde..ffe573dbdf 100644 --- a/act/plotting/timeseriesdisplay.py +++ b/act/plotting/timeseriesdisplay.py @@ -1864,7 +1864,7 @@ def plot_stripes( **kwargs, ): """ - Makes a fill_between plot, based on matplotlib + Makes a climate stripe plot with or without a baseline period specified Parameters ---------- From 7f79e724760a61b924754ac7bc1fda94b44482b7 Mon Sep 17 00:00:00 2001 From: AdamTheisen Date: Wed, 30 Oct 2024 08:09:42 -0500 Subject: [PATCH 10/10] ENH: New image with new cmap --- tests/plotting/baseline/test_plot_stripes.png | Bin 19671 -> 20041 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/tests/plotting/baseline/test_plot_stripes.png b/tests/plotting/baseline/test_plot_stripes.png index 2943429c1f02678a837f5beae6623149d4be017f..da587f3a6c1c6df7c9df18e3dc5c1b21f8aead76 100644 GIT binary patch literal 20041 zcmd74bySsG7e2aC36T;c1QfApQMyx68Y$@z5RmRhNs&@Sln&`ex|^ePBOTH$-FLp8 zbMSn3jNcu1jQh`J3`KD7wbxqjTr-~eJnM~&)FWJMa%>a|h5J}cL=J^II|u*%3kw7O zcA9OLhHsqKqRQ4!&GfD9v@P{elG@hhFU_o98tPoL)w8rRG&6m`$jNx0{@M#`YjZ0u zCMJ{97ciPx8Zc2yo&5oWoHrMHW`#oGYa@TqQuvb%Q7Fm6$0CCA_L0jY4r-USPB7Q{ zEic?H8b5#cs{Pw=OxjqTywB6UiJ!x_miI&TT{1H=QY5;~!U|j-2$wD0YrRn}C!)@? zcHu59IjwK}U%0P_U0vO9%5s+S>IyJ69r>;v{N>PYpwIVU^+&)?n}NY4AuP1(*NuYD zV&Yx9dl#$uBSW34S%t!)v$V8yOB8qiIM>UUFW0`)Wxq(~Ki<`)SINhknVFgEOu&2} z5+bXn=AV~$Us+k1OwA`K2!CmLnb~|?pt!6oAtM7PH#gVU-(RH6cA4PfMKpf?Dz4=a zp*R5n0Tdo3hfMM|YJ6_ax4fK(gOju7$B)aTOd9oFiQ?-kWvhJ-S5U!`kvL{%W<1*i zy2@3qE)^dWKQ15bj5Wbaw6(QgJr21Yr;V?#Z0wN-hql29UgZ5~d^j>oto8D(L zlMQdlSdA`qb$45gx$MN*t&VcI9&%j2e!U}Ii7BPmsJOWJ#n9}`OzpRC=TYM;D|mEt zbZRA*=*In74X@EJ81b4iRcsO}W@`krMRDsk`d_}wz|bJgf8@Ee^wN3zhk^Y{$*lDt zlUgxBjPtMa^b8CWi;HicCQHS^pJKCxueP_a>hD_BWPXZx(M@MW_d`3n_pPrweA5LS{{sJAEfJ8V;YTHn@nVT1Z&TL!M zbxzJ_SkxlQ%AZizJ~FDSWZS?`a&qrc*RDD(CJ3BiE?PX=(k#SmnL0w}r9Zvb$5=y63U*mx(s-Y20Amr=rjzSZ-dHXi`Qh$MQ zzh^)|Ku%5$K~vzB?$1S*)8AFv!`=|x!A4DjMfAR5Qpp!X5uEJzS6NTqx^?S1509E! znXM!ib$y;*^Ep(m?Xrr;5up%PfeHi4cV(pPEo`s}4-L&@4sPxX#>U3)K79D{RpRB# z4Bt0zF2VppSo2>?V3P|>J$zL#N#;c4%_7*%NM5~qRXZ}GuKeuTc{@A1iPcr6zKmxH z{oA|iQ`b2-lm~`}-WHk+`K?biy-i8EyG+FbkA-WQwL#m+DD4xkQ=z$)6%iX7TWx#0 z^x;|^zT;xH#~Hy3X_^rdn;luaJd4dT6xB6Bv}PNzn9 zl(eDgFb;``phNjy<#p_F_`E+8o8(P07j9MCK}Sz7k-VIxUM3fHKUgmE3JMC! zDJ&!zDYbcE&!@yc79Ss94|^gC4#TseqT&;WMZDC6tY-KW2?@#V?WG?#9tK@ySXx>_ ztl<;4U6$=wd_saKG7oF(OSZPQI^RFg3^F}3`u_bpVRrd1zqygJqMY2^i_3!rL{p6c zL}hlXeZ~>6F2yW0n%A#iizG`$c_t<%!l@+Q8Fk8r1B%g~&8c+%{{4C)X&!1G9v;uo z(9pT=x2%M*sKM?T`Z87TJoG1{R#R8^DkvzJ>xlCp z<999BcxNObB^9Wg_xug`SRZ76U}(flmVSJFlx?y#b&-OC?!%Krb8~Y7lfnFpaQ-j) z;8SI7*akSPR_vTZtuK9TU#PD--d|L?@YCgB`#eggu3Y}Fzu6u|@imW~J| zSy}4g;bCO?f3>vqG=jIdhQb1WD=RA-_f^8xwKC`HS0dyC*v&>W3N*naz=fHteU-qX z9<`ssMH47rZ=@(-WIA`jI^GStPWRRuU!yp=)U}oe?bUU4XP+iXyh=}}Lwwvq4^4(< z)zx;sgN(xEIANcA7Ly_$LPNFpH|OU@E4U}8rY>H(kxfcW$t-T&k(5!NS6-o(kd*1nbltDl|a`foGwuub-se zKHVI4k&Mju=g;TcmA|%YQo|i*jf~+($ zGSYnSm)^^fQsRk;iGe~B3b0cRH!CYE@Ok(U#J~>EHZ(LK0hfec8TH`7gO&((Y1^2Q z;bD}RZ~z@Gt(K{&X}*hkJ2A@F OZ4k3=!FAbLv!Ip|@ zXYg#==VqCRCSr5BKpt*?I@F}Y!voe|})`nb^!GnXr-r6M_AL736` z`PVz;9IaQ|8|H_*>%IX2UfJ2%ON9ZQl&a$5;!*Y@BIi(AU0)tMHX}PySXjtm`4j8O zlP6jxCL^&X$%%<*s8XB7B=+`}790%?4PUTi2wcoI3-VE}l93#;vqxu95F@>+s`$a% zdFVC;YQyK?E+}m+lMSALEqHo+8x0Y3@(+(4;i6qEpoiy80{(XrY}#=HyjK26R!Ry1 zA0OWtcQCbn!(_SC%`6t<7n2YVdP>4A=DTi9{QP+nGJ=MNx1F$0? zuUgnLF23UT%q1%}TfIyXewLuuC|8z%<$U`5`MjCgs!BG;Ww|lRp6B1GR}OYoNl~6v zRaJkrwz3wFkZRZq-pS&04Y?yt0(LRD3)@jJ+td*!Y?D=a=hm(Bu&vB5I<7_8DDaQ* z9V|VFa;Q2vE<;^&c6Q#(s_*qmI)PmD`u+QI6xpzFIO-kS1v3uE;JvX>lRIN>{ckYO zox_4TBhNI-9vmwfhj}a#6TlYw8^JPdM(fDGefx&Q=goy)akJDhkVh`&ucD))*A5M-q6iWm1(9gFY=v?-{QA;qrLm%x8lBEmb=ajvp{vw^<_rF0 zvE^WQDhM6LY&y&lW$(hz4(CDyPOPY);46(7oat%;$E73B@_WT z=1iL-OQnvJgRv7M{xYL3H9^d4ADJ|>0;)E<#4)Urk(dBU&s(O-lmbn7;RS5$b0}@_ z_U5h+A9|y%fsfK$Wb+&+MtnJvw=lW7x*}P^EW5QgC6_g1NrHlcjvKVJv?AbOj2Svw zTU%|8b9y=sFK+g!ni=Yb^K;?i2Kf8?UqX3?g@uip-D1+Lq(}L}5{+~}fBsxMJgkO8 z!PYX)2M_u@3P~5bhAP!wnJ(R8&nV(U-(Q=kv+3h~IdIQnq89zfbXvL@p3!~1SFR^V zI~P!sEuSA?%;J-?h!vNV%-6FtwRePJhsT?bRk}ECKqU7pDJf|yN%&@9B3hi}&S5_8 zo;$)KW;J}(WmC4ofr9O2FmHSMtcX<>Ijy{B%KrVDj(e`c(M~IzXaZ`L&d=`P2@DJl z`W+rRSLihcUZF#gMH+}XJC~F3*yptE4JAoZfHS;XdH73wdDWbCK_wQSk|U?jgq*I! zm}Ei-OIAaJvD9vr5%t1Nu_L%0*TBGFVtP7p!%NoX!-o&P;BT;;oSZyj8nP?4!%^2p zT@Or5B&>KKo5}32%+HsH^2wK4*qbGVpZqN;6A21Qtwzfs^&(@)dK5r@3+0=VlA_Yj zq!fGS&YcVJr41R>j#Wr+x{`5wx*I(8!Qz>>frEa=t+~#yU&hN}YOkWh!;SKb85tOe zQINzaWDXI$M(66}NB1cvlC3}@Dty#K*uV7-VSMeSf!wv7oqjyoUhrz@C@s} z^%3s~i3F3!(A;1F9k|l1Nl0dqRu!ukQo~}t97Wx#+JQb;{5Zki>R&+W8qoLhYvAOM^?`=6{nD> z{qDx2;OOYUt}eN#J;*O$iRn@qaF|{U6=sgIm0HioFMGu8v2HE)<1;fe?~a`u`#>o- z1shHkc$L==DrU%~&+aayXRpo1xEZZRt1H#{;Abm2;R^3PU}tZF)mp4iiqM) z6^x2fgn3iE39K&~#cma^;ENrhcaO`)X<|qjI^}V_I>^0ukjlqMHjWZ7R_H!9!T*7Vx zMpv7jzl@#};m()v@dyvabt6nadULVq(nl-E6kN z8sep6_;w)#T)cV}F|N(!AyNu<)B0~AcaU5--)!tSa9iAzBHUhgCkpDsqwJG6zR$tK zVLl4FDypbR0`UhQ2}0nHV>LPh5<09M`>$@z44`Z$a7FAH;X_mGf~&O@hOSGUt{4QEBU zLpjmwwCAl6nqeGB@vB7J)HD!cI2si7^)NXk5z_{6WT4cB?zUw37k10I+pTP@E~)eT zTZ>%F5{8gTC>NVk!%seNczX5|U(Hb7XR%BxFQ%R}4`({dawA+qZ|~PM zQUoC>5p2+$+UKmt8$Sa9(VB1NwqLtuI#TK-oA8iz7jMN9$!WnE)eQ|%+;|k&)9CT~ zH6WMUOPx;AW91G+Mg`cdH>w4C;Pyf050#-!N@l^IUk!K&NI+qt>*;rB=>S-v0aFD9Vta7hWD`I!oe))hT_$uzFPoGkbgfN4#CyT~hDZt&c z8oC8h@dtvNPe@4-Gt+qX?6P8pGMa!RWtjJHNtZGtt&qKOTk9D2Wz@hrpj!QMB~f;) zJ>dO&LIk;*ZN?B7YzJ%zr4W-)keioRoAKretjQC%EfE0J|%JX-mB0C!2zJ0q6h{IV_wCn!-!yWMlH!Dgv zH?5H>wb`0D>&$~$g;O;%g%@_6#4uV*(gnjU+xwyz7_xp5D`wVLS zp}D5p5i1;JPnxGmFJe4g!F?@_l#xO`I8b7hTti`kOLh;{6hzeo4i_ow0a2ZRsXT-_ z8!(nPFu{nICtwF2_4OjjG&lcsW$?=g)7#G#npuGd*}fgjxeQZDAI=&i6!-w6Vqq~4OS7&D)fT-6;?8fFt$R1lI=Ih6ydiXJ6 zl?B`1Ct+}uSi4sLaB&7;vi8%VO8%4%Zz?Kyz4566k=h!nGKkNZP!%xieZPbQYoBTg zYJ!s7T}ji0nAdUBT`twLsE8dK8yhKvMIJuX*_`Xd$d+SG1h`szeX;?z?G6SWWr_=4 zX+QujU_b_I<(1v8hwG0JY#PEk2j9`=QyeNjN>V0tvfo!Y0; zyy*P={1P#IS3Z9H2$&Uz(>DD|#g4|yvC5I>E)}-R!GZfZe2+Y3;ID`aaM+y$98&L=+b>6pc{P+^nsun^;=P*qN4a;bRBP_Z}b& zU>n-IYZGgmn<5q#%-bs?Z#_NFA*celxp&KB$nwE(kN^P{?L!%v;BC%10FCOyS&eZ< z=W=OY4&Rsi4%LXF-?$ljRO!4Wi^5p@H%bRJFQkE7q!3wVwa`W)F@Hc z0Ngy_pC+{=2*$sy)h44i5>rZOrjevn9AhRCxEMb$Xu`CkJd$V_

&-b?1jvDMj0udh0YO0J*%q1x;- zTVgh_4ASm;nVtmBX#Iyl%p2GaKiHP@fVUthI#PwMcRr$e)6yalwGC&L6oEn!NEF*k zXm;Jg#-W8Ag}CN^Gu`ZkpyV>=y*y;m4LRC%A_v>ckRf>G6cv%e`U#Ng1R`d^N>$tm zZVx-AkUTGDX<=C0lv)HFaAt>J&k^I{jfVaVmkEU2oBC~$G+fqm*TDUG0ceX9-^%$f z-i=O$;_yBZ7v31St$QHP%ih~w1i&}5>!E#Tp@}N85wATw8UPPpERSI|yt=AjrD<$h zbX#}U0dPrs@MoSeF);`}iJ;=m&EgR2*C9JsTIC#ukNy|ohBk^XfE${TTsBV*$Y_Lv zL^BwQyaP2&er;?tz)`z9(-uXyzq85&X)01cnp;{DLae;h9LB`HyQj0H_Ri^CJV7JX zR|I~xT^gW9;M$iaCSD#M9-r7u!a3ZI%WG@^(M^CFE<-t22pS`gA%=+pWOqSHNr~Bc z*BlVh&8v3|5J-1-r|N_OfC~gtM!;J@D+qog(sXt_{uR4((KY&dQ#io>-8< zKp*M`q)K;?Ty%~X`>C)gH|rP^lt4%=4i+z|47Lb?_HF9ZWIu3Er-=jD(42{*2%OsiD6 z@*+MbC+A)iT_~g`cC};t7(u=&&_*()W@lGoss!W_GB|BO#{!Pksi47eQP#L-ESBiq zH>l^d*2ZffsrQBsB4i!_aKA-T7YO}O;~bZjiNS~>Juh5I zgpeSvt*upqtH_0&zwI!}AZujzDz3#5!J4sO1_zA1ZCW*u*QeI<~jwifbK zhM~6Z?(Xpiv+&Eg7$0;lWJ)lp2i4RoC-(rb7^rY!2Pf45)k;+vH!=f}O?DbuS|Llz z`zw_Prcf)9n?wM}4Cm^HO=)RqJR~NPy$x=ZBy%VC;K4QmV1Ke=3aw9M*mWr9kc_CH zP5SZUSJ}yTOnbG};R!$p)Ot56K;Lnq4WEWW)p_#{HzvwVJ+lCu0= zc-lkfa%^EA?|2hcGn9ZzIQjYWw419`*iKGn&!4AQOuWXftc=U~8Z!7?KkZ#)hES1d z+>7pqe$z8)@Lr)X#<-lq?Z(y!qv7H>p@Z$XoV@pgmq^myg(L`-*~Z;t57S9ZG=g7f znOwp4s-++%Hu?~ns-=Hh@@)0;kP%kg83OET3U_jWEyEsWMJI~ZZI+?D94xQTpL>$x zg#!GriRl?^bf>5Hj|#u08$DG_dT3#xkTGMu1qW22=}U%0Kv!4VC;d`MiV?P3x4iok z4Gl6R)%Ty@)i-0xEX$VsosqAT^w zD62fljla%nX{BIp7p7BWn))s(^Yz=erV-^l#9VfGUI77@n%k9b@Tz-Nx(I6UBMq*kF) z5UGE;wen2IN_N~d%c@eVQ?E=6PQ|JZZitSdp;07ndadupu1d7(t;MccubBbEz=$|q zPNNT|OFcBije(h5A9HhQ{VzYYq0iO`)d{D35m*(+79RP+GWOfskPvn`xi66~lF(l0 zPBr<|{&?nh>yCJgSNJD=E)M0Q)?8hM-0<$+A&JNo3U0@I+N96XEU9Uda_`9>1*Q-X z5@M5*nyHn3c_{w@`-$FNHT6U`Ql_*x#dP!43ONeNr4DN1+tOJoD;?Bz&6gLti!BtI zXt)N~)-!+1j5-NV%-k+`Eb*amXFxY@Z7o>Gpe-d=*9Tktqi4L8fmtL+T%$KGt!L#v zZ8_bE)>A5eJOyPkp%`bv(JNd6SE%rNXIimb#?kDH%$N!|-L9rpI9+sHtyr_Q3ekD< zM&=5aXUFFU6U#FgTbgcyTeC4(>V-I?(a~}d36J8S(M+#F>Bp%4W!ded0=u1u$S?3p zB40F5YGZa)+S<<$h1{pfJp<*>yqc)+MSbu_eS?&A7Zr(TS~^FH_>W~K;(UZUD=RK> z+eroRI?3Cw3DVfEM3Ou0q}_9I?6m+$i*?}) z7mIkMOnQgX*FVG0Fp0JKMIX_{xZ0B=&AVJS?{#u>jSdgLV#kSptS+-MNRf_oH>$(% zl#0=?DL>i2kal!%k?T2D`A}h~4!fy9Xt28sKEA!`GeLD_ECodpq0G#QOc^l~TH07z zI^h&A^zzo0FY|ubo`XMXh9n%nXP0K+q)qC0?(fwMjcWKMCJWoFPztS9QH^G0$!Hd% zdS9KzisGdhU0IQf_-YV0SV9yB)P~i-D{*`dglUoy1pBoY)>ej}84(ek?QIOea@mj^ z+0a>5F_GXt+FvX+>B|(JYi6n|GgeATNQm3tR`q05PaB$vFeBJ^JHcF)iu$7bl3t$w z#J|HJ`dy$N{JZMduVXf*e;FF)h6lK@X{)SZ5&Tr}gvCrIr{V?+apQJIFwICwQ#ca`Lc|9v3p0Q%Y=#X1sg>Y`7lQ6p5_aqp$T_; zOG;$as{%W0m&Z*^@oWrfn3)eUtKE~pm8es6o25vsRWGWaMk?!VdvxH?vN0%C) z?-V-J|3;xOoinpLNt`rwjNMGgW(q^dY16$ZGqdM5H$@tq5>8B*px~OVt9jMn<^^@vvAwKcTttn$kr+H;pQ3Mn-UQDOhRc zJf4MyXU(iO(kMTCxF#RxUL=IJwLDaC(x6%rY5O@#-Nx$u$E3`QMj=qXsG=&_A9$Ee z9}EYdqrP|0$M9TVKTp8y8WKLjwU^b)_I1ar)6BlbCn&6KE5^Gx@uLJF2B??OPAn`W zn}h_%3METD;oug}y+A6zwkA{i4bHt~9(GpNlYC<(VsS}@ns1NKJ2}zVPI+US4ABDa zrVtTHN&NEVhn&1LKVoq9+HsKgP(2ZseHKfnt*_V%z8qcPyZ});r+WVO1yLjSDl^7~9Mv!a054$i_Xwzr|&{**)`RKmT@LdYCCamkK+{l}w@4 z(I?gm^1U+y%p=OJ^9#9cQRzC5f)*N$q|R2qyFt4-pCA&(C`V8G1j1E%5}H8fM;+*& zHQ-(8GEuTp`5GC=x3!d7R^=UN*2;EtuxNIYRrUdjPufCoz=DELaLBGTu)UNK9}{-U&?s^hfU6xk=6E9%5U;9@V{_%PQi?q@QR zR-2<;da&I@MZ3F4)0wHN@Us0XF$R8?-||R$ZS|`s_6`aKKdLmn1U){iZGY5pJ2os1 zh@0Zpawi10PJdm<%k?NN3GLc_Id+oEsb&89CCs`XfL%emz0pp4*W>j#yOz84FxyLa z%W;8Sg#%sn0yq8->ZR4$gKJ6(LT5fK57uPorO4{(2^wmnRehgCOH7mr4Sk&ZJzh#d z7A-N9_>lH39!iFoRQN~9*Ok)ts1QqyX9VeZl?Rg_p}JQOd1`D49%E%;{o_wZ$>b@Z)XOyB;1kLx2%~&~2>F zmX1jqXy+Ga*3>lI@7GM^Q!v};mgYL@+pd_gyz#I)+n_y-ZKzQ2N491!^BaH9JXQ8x z1ri@rtlwLYyr*||Ik}Y@7}ean>dv}KNp)=|ip!K1S7-lfCP9H8#7hjog*7NEXva

tEJ-n z(l;K;yCj+XIHMQCx3<`om?kOzA%Eq2{1ue^Z2RKQw>CW2P*E0C%CBFs`Gwvwq6&LS zv^1}}T5w!L(Z4Qio7B0jE2>H<$el_2`qfqy75hXyeM4H)rL`r0)zHIqtnZLo61V_4AoB$`aktVJ~V?n(ohUKSzjX56a6AmsIAPJEzkYvC2~1d=uHp23xdQB*r!NQ+ z9+nI|9bt;tIz;!C)ZcIZWsg4cQ#qChsRAy(C&%e$3{^8 znKIqnf2O>S{ol9IcSr{668zpUU$JO?-?CgK2Gh~P>RKx`I45Z@(5G;2tX6scqkAJrye2sk$>gyN%Sy8s7 zKto~XHD*2i`10HI9Z`{y!7(wKCbFYicv}DNoQ8UA? zuLkY$`T3Wi%}91nOCZFyAW`Svcl=6I7j!E&kmecwn%uVZsYL>3?`TJbREaheX@Bi) z%&up<9bJH$Py{Hb2&)xZXW@CX3>c>yeaGqsN5Tfm)5mAA0h$hRjoWDo+NZ-K>v&Ff z`r@WD|9P_I5AiJSIQ$=;U$P8Y)srll|C@Dtd9BmUHO|OJtzi-VV;*=G|22>QwJHL; zkg1C||1s+F*Z(yZH;@1RA3UuDcDhVAdW=}1JM7V~;e;#X%iS~o z@sBext|!<5jB32r>jDQfPk^Zg46={h1M43>PLHy2J+U?LJb^1I7$>6zgf%dvZsVr_ z?Xw#9A#`a*Mn+(hyaM0?2u<$L!ysBw-av4TzaMD1QI+6Wu@gJn%8!Hx+}W|QeWGS? zHUXwQgTlMYgQcmd8MTM>7hn?1K&=A$3=T9sA4ag51b_VK4WL;0{z5vSMb&0iM}76& zz2$#gT`%^g&uR%K{+Lw|7`x{=oX@PyCL^e?=D zhcH)vg_06s83XSMADSA*wNZxArw^|DArB_+&70={-Y;ld1=i!z3=GFOl$7{4cmhsF zOyhC~h>3~8E1(*KsDQ#(WJU!PzI$e)JSYL6SR@p6M*i`( zc3Qp%fV37gaH7s(VzN6%F`Pd=*#8_?7Du9Xo(+kF{?7MEzvS;&&_1x*&YJy@7U*!$ zj6wrGFVf4G$ky0{*4eCiWAW3hmO=)AeX!aNILgC*Vuh?U+#a- zD6%;Gt31BN#T=g^B80)x04PGx<=;y@HO3H%#?f=a!>uhXS{hSNoB=}Ia zc*$(E@us6HLjU=DsjGNAm@dP`7ASWh#vs&T;6<1u^#6GO=OYPa{;DK5f);(T%_0f( zoIQZ604;>Gz~Uq0vc5}5NEmL;6wZD+WP1FWyG+?f3bLJRHdTkenxVznfM9p_>pH`w zHWz>zq+D*#3{YafHFsyr>C<{vWY%0o{sc?5Bw)xMfbW8JA1YU}BGVD0m;I_Zr=vH2 zQ0gGYMu$(mXTPQevr%}yE3>&Ty*c$lYb;=^vP8Ulrsy4{{0k8mcwPA1@KVv%6O7LF~1AP0)}J; zW+#Ne?5MWC{z9REs>EvAz<_lc2=gQyLW&nE_hu-7SdVPe%E}7TRRQDd&M5y3dBZ;` zirxa<9BK_JZiFbf@$ZY?4d2^NBAe7uR!(^uX$&vB07Q*|F@Xn7ryeV#kq`3-$XT_} zdbwxZNB`mDNAj1t+Oh|%&#kOB8m+5DPfvIrxm8FpP-%gOGXb15NaXRs#XQeZ2X2N1 zG&MY$p^F9@jKm+O5gOitegv&7pmtyXO_cD}Zzpj6aX5;) zri5V;5pRFHRVRVDx=n1;D=gyWum4F>^3pT~$10p|K#MQsK~=1P+l$gh z--|bL9Jh=RW(MH7-auPyJ4D)Q&|@fu@7FxtVIGi2Ag&sEbi%3flfVo&{^eq1dg@a% z7L>kucu+kwU#sKVXpeeZGh8B!!wB7aEqi-=E$Dh7#3~KfJ-rGxU_95AJ8YP5FFlLo zaX_H$x#5xrhet<1bm{+j&`tU0JOu3Fa-lL*3eu#8fH8~@f$3>bbivhAqgB+sVEz6y z$$kE91gi})T8h(?@>Z!sGU#8c_%=yw-F=!P8X&HtGMH{~ZEy-akW-K;m-s)O*4+{m&9!eWpW|+V$P+F#q-_DSO?Q<_apv^f`^p z=A^_#ZHUf}>(Ed|G*+m4F8kJ$oR;UHKW4uk>#?cFa)Nt5p$T{YYf zEV}mRbV;mJZXyI`PfwQer{po)!tbuyvKZY3^c2-}%arWudX|DI%+)k}QQ4~E^hAc)p z;LoZ0hO0eGu+#YG5JSfP=MZQ3VNZz`=T1f~c%Pp6_jJ`z75_OY3bHjYiN-TWR7oi* zSg7Fl@qLXfCsI$aAKBWTk$c)3x^524aTKbn>x;6H5y!7zgu(CPviyCfW+YQ3^^;0} zRZZNb|DaP+B1eDc4rP()t??O2!~L07^A8`yOt3E`DqC1^eEm|1=#{qNr|t#t%cv4>$LT67j$l>mh(^+nJqS@ z?Ud!SP)<+FMdaN_V6|`n(s7ztx^LYo$S!k~=$0~dvN6#G$g!N56naNA?{q{($~YDlXC8P2pnl^GjOYOgDFNx*nCXxXjCOQK7U?k-!L?1Mqp z`HHgw2gj+)6>{<)Vn>1&Oe`&}b;I>o4zLBvsxDui0GgV*)P`$sxf1Tyw~hl1qHc*h z(99u5Z3af=P9^R6GuzTLX=JFZEGX11vh_5|;1Lofewuk8^?+eHIV^<2@lH0RZys)uq6AzHi!X6o zo;&~Tn@m%acz15Dn3!bfU?e|N!{pSv;E%T(U~bB4l!S6m%vSWdLGyTF>RSk4oV1kRzR5wei^g{I1qY+16yVq_SFT(Mq*JoqVYq#pQYrhk z{8c3y5(WMI}5dT2oTu zLFTiheMh`2*iDJz=R=P&Tf{}q5fT!`@xeNv9df=u>v5BWm2WS)4hqlqF`;79FZ%#7g-7eJr1YTsWUS!TB)8$0j;DUp~Fs5CY^pvs?ef^`=@e z%EmgMrNYjvwa1KiGKaE}SnieQZJ-5;u48sU>CV7tn# z8+1(BM!ez63Z|M&it=CBGrMo!!Gv)YwwQUZoDbdT`SC-prNiO#CtuGNiIk^uU*@Np z-pTBlZOwOcm2}T6FFm^{Uun=2I{UiDrhnvbMqJjCoBiWAu6U+n4_^zR+35AY zTw@yjKAqqE+(0SnY$#`E50)j1TUeizDc?6)teREMn9;~il{4Rx9iN@Sc zc!B7muy_bmhTZ-BNku(vGylv?21x2VK+u`AG|7)~VbkpsblsOmD=0ul16_sIv5_mk zqtSlEC1PS~PR`4AV=nBVtP3cVeEe7t6irK9>gwvxBH7RU8(ZS_tCv<+P=8lo7f(|b zBa(kppNp(#={pt;kM_4f>&qd#clSLv4KME%p!xulg$`&S2y+`@|7DLrYI@J0jdF8! zj0$07-MwTp+6Zhygu$l^k|B^6B9g{?&;KHD-s)9KNpTi=VgrdcLWSkBUBUv2B@Q_Y z7NRZGC_x}zrSo1wuMql@X+psa?u*MM7jSSqh2lD#B!SvuOVJW#*fZ;* z1QbJHHWC4Y3dl&T2+14td-dO&=6aI~X0l3V88|tAj2~oXt(NtLv&Zhu;-NqUwFZ2p zc+DzTpzU~+m2nFR2@wzz-zp2}i!9R@m#Ogd@;Z;Y&c(GKXj)WUth~9{FKrEdOhnFr zNZKYrBS!>0oJkZS6Bpcg~lq|+HM zS__7Ph|IEG%*}HUEvBaHUgAKNTNNVYLR~{FFv`}fV(TLmVwaHW1M;>@C`7+Jx)0KB zQ7AG+K;$e6^5%g8<1qWB+$InPJ8ms-t@=>6^4hJhuh;hWDu66r40LY_h(VNZ^{FC? zSAQygjp0g{p*h=6Wl2ggJ0lO;d2Rr8lN#Wk+hDHsz?5V}*ipkHBcfoVoy+-s!v&t@ z3Pn&3cjoFoLKAR?VmEKdFw=o3UhgO4wQB;_i+z%izKP6+1_xg&yE2nC-+RXn$aa?y zLeCX$+e<(pO-1xCz`;&OM5F*%Jt&Rp3oo0AE}ERwsn{8WQnLR2EpZTvPKRldk4PD} zMi8{}+X-PI1XXN+zkx|3muoi0w+=XxALvrBKyhNRG{BrNWY}6_oTVJw-#xZ<8^M&( zQ5T7c69a(I^oGx6-xoA=;>jimy?G48o-Zq~H-Qr?tD=%d+MfYxE>VgDkQ>71!ru$| zxOjLZ$b_e-r!fvT!8rQ%Z-W-61a%G6c5y^8ApS)dULeu484U{$cUl2_MFcbgqrbwZ zAQ7fc7}8z|*KirN^9Ho`0!jx|#1DrGU<1FD!>|${wtj=qI@{14kp(i`u!Foi9;^g+ z5YMWAb0MS*5uYI~Uf4ba)UurPys^?!B-Pzgl61ddA8U=K(> zpk4wxAv^49WZ8_L@G!-E{Tm{LxR~Egyz}cQ$;sCU z&l^auc|$+)^e%t}A;d)^-M>lGY=OW)?&*;?8E(m#Fto9>d<~i~jl-<6l`|-{3P*)x z>Fw?9Pn=dsQbP_KGsrMN%Tlg5_@bDO$U@*L&D8@7-@mq22t;D(Z3owrOQzuh13ZvL zJ2o2!9)Wz(-_)n%z-i3R8-NMJFy@(<_hmETMIcaUiQ%V$2SyLrhsx~HUh+@e7nLQ( zcV|{sn@51bwFaP^?o7+)IiRan?93u>`^YZLs8Q|<%>kaZ>az&2-?@DF^cMz8)i6+} zYCzwJhlj`2Tx8sj2Pz^!NSRBc+~iimFObXye29$n*Qjuu04rz&f^e=*1GZX~>yfB+nI zoBf}Qi|IJAQ^z;Y!3#u<=*eAQ}h z0fDfnbF^w+J1%|B>iY<8a=Hyx4S*vK+sgrWh=W)5_Vy5JHv%()!k=foo(Olaz3)@m zOzp1TmG`W=jW_@w2y=tyeg)7}-N?ulWY41Ajvb=5L6TmZ9OIf61vb!Fb#i0|AWcnO zT|(dTb3HvM1fMZwy@3`FtF4kzfkFFKXOMQ&fFGE-b&S>@+ zUcAjr7Uv8DBO@O4o$4Wi0hZwdDnX1@st&FV43_rc0FZ3s^X*=9+*t_#ktrhWMP3T+ zD;Pv}wexvPRDDlLw_2}it`0h|HG=>NCjp~Uu5e^USkK5CAllniZ0E8<@R`m`0%;gR zIR+dZARD)eZOH^yh-J|w2jJhv$H(0PF`0sypxogW1H^>Z*4I;a_zL@mSxL)ZkBT9Y zIq&&-U}|b1@H0TWF+k?D3*K@asF-J846hv_-VD$vgfELo9uPOR#Ac+k zk~W+=TvA#JQIH%dsNs;mi;4<}jlBlvqk^$4u@^%1@GnITjX>ySP1RshW!MR7i5a(2 zbR1e@qE;&p7(p&z2#}^AAW>i2`RWjAHo{w9?8~$nVMVp5+a$F`C2_z_8nvJgu>d0N zB#0nuAb$b;&c`w6xnX9Bv8hA4A;J`wlx%<_%2?(qtYQp?|G5mVNWf25%&dj$VZ1;J zg2=op9Jhqb&F_KcZn4xDrDC@d{dOn01iTv9XZ=qK8oI0Uh1)4*L`@U0GQ<0Pd1ev(g{# z1*z%nB?jcW8Wc)ss3`+UWEBz5?iva%MS?x1XJ zf@@Q9+)iB8O0CaR0iQf-CzEv+jAa>dLZCrGn5>}XzYdp_z^OyIL!^Ws$(=p@sWUmv z1U8!0O<3xLg0Au!<*VDYL@_>{mHK)Nyz)4JRR(PVj|DBgfP{p^zvE^m*(~WW>*~jXEe;U9(YL5`!Ki${EM)Y1BTwIeNWXag5EZ7D^q?{bZ za&T~fwgKix232n*4F^#QCK{T~XkW=hx07SvprHP@#yXIk8h0%>G&btgdY*>_raqe2 zY15QGZoV5jy`W0X+E5K?+@^t+4x+UI4<`e*1Fx9`xe*!)h`)t`iy5h$vnz0y0OUQM z1ZObuC1bQjyHkx<{EyVcv1*3n9A+{R5esNdk=r2z;YlqP zRP;r%$2!cc{QMR|AfyGiE^sX)2M6vD0`RQQ&?%GGuJwF&EciOE0Jv@_|47DzuC2*QJ8}b3o)SscLVHdWFfhBs~T{(^%vw};HQ z)^RD=VnG8y#^-GBaS$tWwdJ1qzw9wT2L8)cw_J)ox8Q&4n)!!yw(S1y)K}wg*%PXchE>|1z}!I>2#=A3%5@jHNCZ^JUNBOI zY7rjDMSY8ynA3zo8c2wsZfh(K@UxN;ayJ2oPJ0fe+8QYSeJ|2bIg_ zMHp?ZOPuL4q_XYmR;{Ygbm%OzGmFv#*&~`*6VJVS9douSPQQaUq{$N%JEL`w6aasT z{SqmuE?5a7;k^esiD9Uu9v0ONA62*|1G)~>D8XW+jRH-)yO*#tG%Cm;yUD8(G+Qn) zZgc2XE;79WDq*>7C%gRZarHBdvyb3XLsSJ=eNRX>^a+|JpGm?S| z0p$x1D{4zrL<9>36@-xOkfEWW4x~enkvtn17=Vt&LthAd3kBUTn>*4SV{n)Y=}{L+ zNRqmbp{s`Et)y44UakH99SCIi1YQmo%bHb$AF*s=V_|g|I>7|rfN5gEwF`)Z0MSQ4 zoYD&likyk&)@{Q2`>H6DdvvB!_!jXX4N zRkv!7pax;*8FxfNtT$BHSAM?p&?@ZHr@L@5$Xriq2y(p*P{iZGrA@+PA%avpZ%Aw9 z_;_ao0usjQcVKO--}fS|wBVbXQauXR`EkMv-sRIx8x- zXENbiN~(pg)0HD=y$ORr4awFdKHv8;t0wU*F)sO3_VJCw?{9BQj6Z)W`uMiw#?w8w Z6ny_{3gU!k;aWG;V^JxQEFrDe{}&>4HY@-D literal 19671 zcmeHvWmHyO*X{-s5R{NoT2M(r6ai^z1wl$mS`=xd8>CS|Bn1Hhk&=||6bU6ny1Tpc z%>DR$@Av&VUq0u&B9GM%H#uh|qhW&}g!IO^(-Rc`A@7s;v@BV3Ej#MWx4CMNLQPzO_dZ&7xccoqf zB|fVp3xn!C($@VGb$g@!l>SRUj>|c6p2pw3J=^B=>~X?yKh1?JJ_s>!<-55#=t_@Ra8hdi!HPh4ZAuzJack5 zlaiC!x@nO&Z{LQ-#^RHZkVMDCOniS!70qYOP<617iMNV{T)cR3Y;G>#inyQN_QIct zfyj*F54WEiOf&>`^z=wsS>2@MF$-u4p<@&lrj-t+zGrV=)OBZ#^;ZNHuLXj9jEmF# zqkqK4u3bi&EFIz(5kbc9u<{1`0_E7RIDx)`7bGujY&y#wERWWk=sG$&9t~KuM&6kF zol3a1wRKak27wV73>KM}IB&Ba>@2??^*Gf|Ki2g)_SLC$z9%k@6&xJ=qoqYmUY=lL zVq$D!BFWyYHHx#z<7AW0*wi%HVRa<;+qZfcvDH`&MoDSu_*`f5O89!kY(7ECh@a7HaxpU{laIZYh(u6w&k}y3s zp%cR;rlW%kVc}GqwKA$;9epArBh!>qgZUg*KDHc4htLLsnFKspf9UZ}jNa^_VhYV(+uC-N`U8h1lWjS3rRYFhi(@kz}QoDt2V(l_p;-+94 zgvR}>=snUCngw7{$yPH`0-3-OCgT^2-GU1&IR6Mvy64)%6iEgM( zieq*#oe0M(Z*SkcygcL+v(MD5ly&1~uK6GjYMc5@$a{9se{5(tB&!3W6nuR`vg@tonDX$-dYDO+CnQ^%vIgT$Zu%eTL-T+Lb@7wc` z2ZIxv#9(V;;^Mxcp`mHo0x{?r(}E&Z?wL-eESTtKJHJ)+NkzbUD3p=!sFB( zL)3PmJH5>KL=GLbt$4X7Ap_pzRpYECHw)zkmbbY>$=^4%aLMeR^zZNrYZdB~J z^1Xya${sv>trmx_rlkc>$a%}NwDd0cH|vulr>{DdlyvTUiAl@j_5QG|KHyN26%vH4 z?Crn%;!}*Ru7-d^jvA1*x96#-sL*@*)Mt0Zg%B*Tuf&?MuD(9>&}`D6r7{r)M=Rc89Tnb{Pz087`?TzTA2qA>JmIo0zEyEp$aDw zKR>^q5m3(I=G>sq7vpk?(j52E)o<_4A|>VJY$qNA z1uwvxdVctDdBl0~hPt}?*wOyx-b}1zuFb4GwSfJ3DJiK3hwevvT$or`p%D?j@j@<1 zxw%BJ1kzwBxzF2a(R)7s_3@sJ3@(^}ugJ+B7M7z2pb`U>rh$dZqInaKty=4s6e0W7qZe zQq=nsHZ{Jc$g^8~jNXKTVN#eO;1k`2n?*H^EmR#&CPw;6crcuHkd|Wq9uaOxHp>y9O&CO zZ_b+yl@M!}+c!=$zFr4oPBU{j*m}_xCqO~NZ(9QOv5Cd3n}F?mKtJ$HqLP zqc2hO+u~L_Z!1x%J$m%6Ya0uJAm9}hMFFnPJB&#oNj8dulRJCjI8ET>-8&+=1Yzw? zzkrkJi~|}>h729eLX&H@c6N<%0uCxY!<8;->!+}#$T&DSFg;-jQkBYK%iWs}dv@2U zLbuxeLBccG2+GmX(YQ1;G>{HlkbH8IPK2zta#k}hoA%}EV|w1ab?c(8+jc}3|9+|< zJ~{iNpjHvGQ?T9eiL++n3x#wDS6!rAG_1d2p%}c6_?QQS&1S7A6 zAf}NJ(g~|H5Yq5}<`WSV6bu)33keCK;OFPB8y?nvp$p+*wdzp*#W0#AAh;A16cGv& ztYX6MH*VapUK_h7c{Kk!m4TjKl&k0c?Ea5X#kE zX8lTY^Yeia!6gM@5c8CgvNAG65SNqm>-|_4r);=f_pDL|<2bd;u0ozGku3F=)pp6r z$qDt32+_#Mh(Uk8(FE9*QT*m~99yGJq3ygAB-_rnl< zLwP4d%q?nWeX=P@)I!s!E`CBnLJZN5!fu+&MGs?n$dCq+(<9}Dlh(sC$4zt| zDi!t-S5ILe(RWNE8BY42C#3(38}g@fXGQ|9$wr3rXKx7KYj15OL3|Pu69={-MI}Rg ze0?R;-EeCaIs#8_sDC0+OewG!)|^{hWOMVEj{_$Gv0-#n*J${6s)BE9tXBHjGMKI| z5%%e+Q{0$aMnMeMPo-1Kv=||0WJYjtaBw14A3)-YhZxLtBnGKVN%?hjbzL+cEK+rm zL!&3s@I1loOP7b!S!WH8mh!CPTW6giD~;G0E;84$K0erfSZZ?1*xK4^_@_hS zFygkW9bsqM^Oa>+kTafSzceHJ5%I)+?70PtNwwMJ5aEl6gn>bl618UsxWnC*Vdehs zO-<6sqjLE;dLNIY2Lsv&rjySuv^%y*>(t&X>-7KmDc1*E%0s(y-Muk$6Efb0^@ z4v|N$S!0qqymAj8`g(hNcTq`2tj~2aqG43u*f?#caPKu|6%o?VovuvAweRZa;*wdw z$<9tXFffn=Ystw?TdE*(tb!0-xpIYd@rkx}D|@zp{Squ4WJVEXryF=%P;n(@W!-=m zRH?GKWjP|)=(JsZdaOVkRo$tmX?-K16{?L$SY^B|h<3MLFG8+`AqvGq36s181gUihx#^XK`P>jxdB^wGCp$VjYoR>) zlMPeTFkJ2s`r$+RR%R>YcqW!ES-RE26O)s&LLTsxpn?MKzP>)@tVKwJT9$J?&D;rS z1q0v}_Mnv6Swq6zIHC%(_ex#C6(jjtIe+3>!Qjt+a^|2DawdfZZh#?%#>WT4ETwEW zH#Zx1eDYF)IC}v+vEqFFuV3<9*mjGuCOuy-L(<`0eR@KLLqJ^vh2UElK0QA_WoKun zdXd@Xo#jEuq?dwBMyuRX23pfw6CxwYAzelNYy)JcBO@oiI=qm{tPeYkCcgOnWp1=u z1cGY^?ju>*N zU?Tj#eY@>`FwfK$a|a*NsaHisJZQ@1dbnI{GAhR^KjQGB>6y8CR(Ezu{)TG*lja+c zm%+P;!!K(n=p&k<$ySRCOeaSwU5tWrl7vr>Y`9-^oq@nv1L>@yx_SUS*#O)b)V$Qc z5=02`@n7@8-aL5Mq`T7>nD*?E+hIaL8nrhL;diK;-+lNHNF!jcx3@M9?e1it{mdM8xZ+3@9zp_0@M?MkjqYVrO1+v zWM{m5L<^N6#GddEbXU$SEidaX@_ptPdAV@{X`cFwAe0>s@_v54ANDCT^EEj3Nm%N- z3a8DR<^$YvA8tv`<8+P`z;fxYPc$N)$A`OW&W9AoD)~uw-9vjXK-6^w)C7Go7@Rl6 zAQH%A>LPDwp()VIS+V;Y!qbK3}f z#5aZvR-EadP%1O!fBVL^{j7e`YU&3uln!#&PReFelLUo?Fg_{GZfr01c7DE35}~{9 zzO(dUU~%u^`%wnY>fuX{$Gsz9rK$#ypQ*O_3N~?Wlqr7n^*xKGoNuXk(M$mli+V`K zTZg_t76QSj3C8(nB$C79wy3-G@5NiDJ?uxj)u+DYN;sq50;cc-HK*10May2N^q*PSGC+O&2LVC}({TokoCu9s< zC)=?h3m|WRrJ$P!aJ4N~0!d$ZF>rL&kT#$3Goi8U5rCP7B~=H0~>!_GGhGXfeK zB+zVhxePKuVlFN&z{$inXMRTXXqUcB>0bShRZlgCt;1B;!^0qJl(~PucDgl+n1sX|_6wW4 zXcpLeZeCuk-vz2q$;kkAQHlbx3b}rEmw6AopNWU)0}JpSH+gulk*80ej(uwvR{q-h z;|GF3UHHVhe8uM+Nk^I@Ba%-ywK!bP%WuDg%gf7)#>?BcAxeL&0f<4**f;>-0F7z? z(5Zy*NOn@smRUUjb%GSgb9Ac-M0d{H?-NGd0UQ?UMrOkqqw^OobXK{!pd}iN=@mea zkYaWuOA{8rn?nIu3qf^#e>R~WFg?hmt!!;2baZGhl8~Uq2%6a)hcBq&5fFS&O>htT z@k0hJJJzXf0gw9%ILxxs@t5O+#hkgJ(%XK7bitCZ$rAyYfp`4AJXnlECQm9Hi2)^1 zb6vP=F~kS80A7|>u@nFV=p=#dO~BjSgM>8?JYz3U4?vw20GaH~C5Mdt1wlj=6bPZ> zg`$C?tE&qIr&ZoWI%HRT1(@=aT8_@#pG=A;Po7w<4BcHCDkXuspVji0%8Xs|r%wn5 z28NZDl_;d-0OW=E6Ve&1j#QSo?%zVu1OQmV1**UlC&6meD;&81lZ8|v2h0SbM4sci zL7~Xj9_v?)Z`iO0{2&kFw_EUmH%u?Ewz8^$jEAfhaLmSy=~fiD-rL{5!p?pkp5Iqy z$Cj;A8C2zVXi(#c!FtG2MYvBDG}#oQHl)gB(14pl#~87@E6)%L3 zs27?<`J|KgRtRT8HWeCq+2>N$D$?0=7m7 z2@x|>v6fXVyFZB;GM?$IhBS~UKHRJ)6bUFEb=TtK=MRD;4ef7?udHy)w0%h4v3`I8 z>@dl70FhjxqzpnqPe}1u;%~U@ENP_dAU-)cIWmC#j*pKg!D}W#9g?J)`S@S#caiJ) zw2CU%L;EoN(n%CU2W|>i1R@ZJrt+!4ncUo2p}Bk!NCFAIH@?^8JqyrEvQK* z-G$xYMX4&UM@lK&1;x4>!2l&8)c68Qj8>&krzQcIY|!&n!$eoNgY#Q( z!HH)-V*;V@eLSrg{c?I+dBeHC%B>uQtkFgQiivvuF-k{PS{0Dy z+({(%Q<_OaBIg5OQ1bDDXhS7M3+E0u_edXOa{wzFlshLUCxGxUZwBLLR6Wkod5gld zF#AdepWg=r%sosya9ImT`|i9phVk;{OO=UCwVdyezd=rT)HS*YjYS5A3>Gv5Ih%t` zb!4crp-5z7<5@^LP)IG9`YtBy6sT@;GW=yE;B)mdRIPJ`c?AUx5JX=C7%@3k>y4I5 z(9MIG;RlIx)%)4ui}zKbDE6(&&^I%C1COlQStlYU4&7Q8n1j3I)P}Pt--M#>LCb&(>ctFNZ1h{i_Zv8?k6x{-0d%(jPT;GxKHXEoB3CK)Z(KO!g(m*j)}!6 zfQY)U+<3yjb7N@Wjh8zO6;k7U*j%!+%fNSs?yL|N73owe{Yf2RuBG5bXY6eZgm6={lkF8x&+fCx zloXr~{&|v2`nq~Axurx~(=A}YuPC3#skwxU)0 zR_HzzCH~x8{>7m~9enFs6>26zxH9~GIt;&~VUy%icSPgH>hO0FrYgzXISTrZdI|PRi>=oC1i5#p z@xA%&gLtA8eo4!YuHWR2G>?k1!Jj4gCF$RBKh4x4#N|UEUzH)htFvTH-_ZL0euWd> z(V?7Ipf6+o>UH{s5Zp6c3D_j~{dW`M;yvA{m>6><%vQfQ&mH{Yj@IO3U=zY9?^O_6 ztC_!9_-EJiv9HE=J@cr{j+3upy+3zz1KbmDkx4y8<4LPvB zStR>;Kd>~cJ+`cLsmp^>?zMTb<>b>dzCsDDLR8nU$)#$PxA=4yVIfz^su8GC?4HskX&)wZ}ihTqBqb+o=Jmwq zitLOxXtf@vkID~6i($+Ar$;1DtYX@2YUMx7(a+~<*{z&Uz3=F-MXKHt)wuHc>ewtD z!&+_Mx!k-$>5b*9>>4-mC!1UmEv%V3-)DwaVp|u8TV_ckIL_Oe278UsNP{SoPKk zpIm>dN%vEgIo9VP`8q=U=LUPH5B&SH65VG-NUF{%M&oo^S7R=GFD#hkBy zIVmF{u0N<7$)H8(xt2M;oNVD#QQx_=!}C~1?tZtfkfU&r>RGdnZh; zJoe|8g5qpO+0{P>G)J8tu1h~KuhJ>GMDqOzXNPAyH2mW$jo0oC9wH~W7lIPX@hJ=^ za8_%psv^wfSYAZuJteE67CN)DJ#DDz9Y}ZHyYj7>Twrs0blwzGW6AWO+jhNX#PanJ zSDc;g8&A!Zesh`kJ4;lMY`DtzNKdpr;_Qp1Y94qxGm%Ln%{7urw)FObk|9*c1 zu93{MOK;9wv+)wz1jjcqWl`4sdZcU5?nJm|zNA`VY3hC8T=;HdjW}{x{&6gmT^`rl zI6c!xRK=V4xpNy2U5DDv7GdWLa9hq&=#4in?v*$N_W3UERZ9`%rr!{7s49=~B-EzX z%gn(+WZc(?Ugt2SebhCK;je#uy6;YYY$M)T&0xR(B`GfHZOrm_9tW?iL}80rp_*6t zz3Zd)5A~m`a5L?0tZz4abjL@9xxlQq;Ras58(p+bbUg$n!Y)tn#;E80RYLX_qfS%J zYEiy!Hs|BASl{&O)yCPAtcdJX4aMgw-?zMMu==_=V!~eetG{GC$H<`HJ@{pP61zUjKo+8PZm;i#AK(DgIO*bMEP$PujZ8QSxRl?Cma8JdxMn*^&hDW&W) zecZSY2YhSOUl1WRiT8;poz9CghG1A3GraCEJzLXQN%v&vu=dvt1ks!9x-9?mu;zQp zZKP(PlWynjsGr+pvHtQLA(ETBc!ny)=&b%8VWLuHIeO?L(s`|>Nj z7Zhid*%UO|x=6@*ylxgg(#bZK%QG84hhI*w-gmaPyGuD;#7XyA4~_l{y7RJ!deLNW zo8kJe)UWha%|2dh(r*7kV)ju?f08uFnP=od-uu=@!Vl=**>Y;qVIP_lO~FU0dJ-xV!!A@tG8G{}?38XuW=`l`4{ zRIqb7Kq~9?jY9$nQ&^gAeu?F(xR~bbxa-SxW8GiBnl|dr^HFl2_cYCPw=(}jc!3g| z*_>%q%KhtuUj1aGq;LVJ`0;VaCq>2b>n%1(7IF*OE>YyQ#v|)AB?a$ploy?slFjil zp1(t6acxHa#mIv;qxz)Y=!u-$Z#Dftsjiork-d5xhn=LF>Se6LfD^Q@_FF1GY{KaU z7>RjM<@e^O*G(~wJ~j20^nyp&NrC<^<*J6d8PAv><0NIBigMq8q~I$xzNFTvRlzo) zv~tL!lxTxebLN44|KvCizL=tj#+c#@lI6d~jUN{bebr_bEcdDZCd(I7%x3bNgESsv zeC<3lFHdK?mk9|`G12A{+(kNN8cKeEt)<5Dl9_Ie*2Ls4BX{ToRKPR1UV5N znfon7>Q?5T^etp*n;|V32m8$775IUCG98@3-Q1Oku4XDLGUs?*Xz}J|0uS!y3-9bK zs#XKqpTskxwAw;X&5zvNIx?ksw-!1toK)raUgjJ{e75FrICcB)jP=cJjpo*V2^dxD~Ik>)4FxHJw{#epQO&R(NCART^t(Q`v0}cHJhW!K&}<`o=uE zEZ%vAo|*=;uCt;NrkH*7%bq=r>noyebL3Uo1<%xwB}XiNr$mKfJ=^)KhXR?Vv0rdE zE&Xqv+RrU?s`Px%!Z;|jm8_1;#)OIG|02;i}CYE zYUIB&y*aSI8>?wOe&g74N93`^^GGlqedfX62mjGvB}AS*6B`46q|Xsq>Hh8 z)>xdCe>M$J4=-nR4v$a@_vAY}rG?NEY>%7Xa0yEZ^W7-@-CReS$PBbrncLXkD7S$D@u9Aq8-?LoQlQNcBR^PhlkD{$CT68 z{E=UduK4+o*u>oHktx1oZ0k(Jw_p>30TZe^-6mY~LLuXi#3>WQ0&1+1{2WXzhAwN;M@MK znHbLLdCC9%=2Z+R(*Az+s^|Z|>Ho!Z>S?b^=*c<@h!((8Qupp*K(n+nS6}?n<;$wu zxqnSy?pmH`4z~P9KCcwnCTkAXm;nOYu{B`Z)&&hx3%&xByrf<^8oEm4|tv zQMuoqw;KS?%Lm;IVAi|=6cD950#nB`KAtu|Ki_I|=HW*^YfON28=!5q)<~f{7N#KV zm7Q&&N|^BX(X3u)%X|tgM`-~#+XFzl69CIwgEO0+fQ>qCd!rmt z28N{GAb8}tD=6s?Fx$LtwB1kx$jZoA+CL8nn6uOK;?RdaK46+ykILA4UTFGZ|S~w&GzR z4m|?+j`!xj``KAr{gZ4oG?;iSN1j-e|BU6WA09s5Cj9qm9V(b99_Ztm_2ptfr+}^S zXS>Mh8|XsoL7(-1jPQlpHj-uwsCgaC{!O^x2t~U4zcE2W3Il;o(B5KBb*{^5JlG0EnQUq=mq< zWetGS;mby{-;+F!KD5Rv>~p% z&fbxEZl7s(!Pvy44oG5v5jlDPT@qN6K)3|RT8~zC+n}?9;{BTx!+$+p#uf6pKMg>K zz6qRKzz35(PET@yE@WcD>VA6c01Zq%;G_1HJ8%L@g$-9%Z2l5Z6?*#5U>Gvaz||5; zN*7~z%!2^gXLWl|OHAwqm=Yr!8$O`G_W2c!vU*7o7n`h|Gdj>Ci>5fTwi}?yA zCFNeb$B8#URhS6MEP&NeX(akPy18q)qB+-r+xZCkTK2;Ln*lfX0&?^AZEvE9hXil` zuBi?NjKT+2{MoOM>C6q&u*+^5em@UHMl8f)xa=$&TPEwJD}QaP^1M4vWy~HdB@j5z z0MP&_wU?oxQ(`kqKtRps$)pgk2Q7NU6Cg}mXcq)Y%l_Rv5sjFYF>5cEJ4&{a{Lhv} zEy{X-!w3je>DtdtOoH0l+IR*34vE1z0>(rFko1qAKUK>)?%ur&P_w1u-(g5>v4WGl z4kQhgDi@Ao=WP?{gH-O{`}_KxZ1uBnE!sXQF6ISi>kr*ObQ_@LM=*2fe)ym)obCUt z_8Q5S3D^J@@PdD|w$=k^k9Xw^tZE?4i9Ue;(0SKTB@F#{*)c`#$xP=D>;@n zIMv5HT4cfD|BjIy?)a^^SncQ(?eAr(akeFh=n?}@@rtXP8`@AnTQuxl|87Pu;b&Pf zphai^Dr$j&PDt!6SRGb3Ioc4A}K%bGXHnzD#z$f0TMSCSR%w9 zaTbU?_Gkxo4kn)L?VEqca~3s!^%AQq(6M-h_R!?x_+h^{fETqoJ#mBZVF1(_RUwIg zE%F@dSWHb#QK}Qb;8M1>+(1X@1jrkD$cDgB2L3h#)Bh|93{V1EcEC2QgWhe1GxOi! zZ3O{>eGyjgw%JHUY27k#6u6$XU22OLI=mAc^zYbkz%%!NrJtu2yaY_Gn%M-8U;wfw z!JP!@nsG8S*Ptw4gGM~qm^3>(J9}XN&|MB;rM>&t;8e*it)p^*jbS+TBUC+C599Gy zb+5ZNGiWOo{0y|P@9FDbK}2!A;H{xydtnZ4PZkpM@AJ>#+Xj9sD*6uZ`p*uBmM@?C z5ucKh61X1(nn0+)1sHU>g@suQi5LHQW~Z#4m(D|I64gNO4M>&$y_YYb4T8U)|DR(1 zKb>KBWG`TJ#HxY~ApxheX$+BeOwbQmUg8*fSoLbD{2by*MPYx1aT^bzXA&oWU0s0l z;e}dg><3t#+{>qcH$Pio5a9I)IU56YDWaT*%2J$|KPSiq*dgG-92BI{XY{j|A zkI;*#k;_)wU2?k13oVfqCM@KZatOnz4kq4)@>9cy0>L=DM#q9}3O6WE&k6pJZK%g- zL0crs?5PQkRhCX)(*S?^klwU57H`HeW!NGk;_A9KQ!8n4axi|S;EHMreWWzrGerYH? zA*z>iNlYqZ=}AOFyY+L%(9n&ekJaRC<=63n#O+X-YPYi**CH|K*@SeKt1NZG zwFQuOv5Gx~7gKDru4)haFWFA>Y%Bv{ZFn|F+g~lnJW92Dg>L$;BJXv(-sA1f@^_&x zr}r!327Sv}ByE%oR$5|;GoQb1T#39Dl@}_b5ZRt`;Ff8c{1HoIX=60B#O5$);*NaP zc|nii!&%#_D?{%aZ^=SqC`w+MAxnbF>FW}YK&bs<{s7F)7xh?fjU_Qf_GIJziS#&` zhwlM9=LjceGZnBiEOXU#iU*9nnP>LT(&1L8%urC2xdW4A(Qv5uA(00oL36RutX_8_ zV?a5JsmEf`QlWl(-tUqaEBrtf2LvfMCSKXWU4Br;1uru}`ReCxoq zq0quWen;G1%OUoG;-$jzxmve?YU>LulE(>r6*z&sMzgLa`AZ)~?4A<7pD@bfc3N8F zxZ0>IshPjTcOqOSJKbjxBFS*~9=R+haTsTU0C8Zrghp(6dj6!ai{QoPKsD0Vg@Lr# zQ5z^j;+NEYiom_Iq}u*W=Mah;w;V`sNzrJ8emA&a}; zb+kFWur$)g)ZrTja=(8U1My#dIyQ6cdR64z@!&+a_n{1{D|2jXlW%@s{h^HW5$+3x zr>#sY(+um0Umu!`Gxk?aT8s}2Q>VDC3U3l-k5qa4*lL`*0Kxmu)&9MUzgJoeqz!|9 z#T6t?@_CuA=R;dLGEQ6(hmb;m^L>QDwW!y{{SJH4T_36W9MAMxGAER~gv8esI`)@d zqw^O0qAvV$!znh2tGP+L$KEA4q4tkNzthL)m~3}>j^TzsNk>&fDh|=4U}y6_fbxMuL>)UtWXz`8H#Z zOt39RT5%A`CA|ghp^~UnQh|@zt}E2??P9K->hI}(956e%(>zcfF0Sy1=&N*%FiPRN zmO3eAL2E&=eL3Qt-(vsnQyZo3$HPvgvlADdJp3tY;KyDy>thoCH0MQeLiLxp%2%1q z`F7_VzCq_w+-dY$lkHVT{Ple>Z~E9Ul3i6BqRp#>#qmO71K)tz{JVX_T2|b?a+o+I zqA1LgaTVV)NyO=`sH%cE?v)H@qOF;)TEzrV2W2ei|I1(+pU9 zUceMnYPxAV_Xq?~C=&^EQ#C+KWOF^A66(AHd5J#|acY1}4x9y$KS}^0gX{S(YC(4A z&r%xE<-N&Cy~98jN7E$`;ev{f5VSLnr^maY(b4J_ba(DJ{JCP;NJdIp17vNq5(1G* z%X0FUnG0#OrJI|Z-avF|Nf;X&dwHn0q7%~c!~x2MXhEl2kmEkee0&4g$+bX>1VW(E zj-s*+({{zRt5?%1M_8SQ|8!SibhXnUsCEa$CU;`)m|}w%4&{?wy?O>@No>U-F&_n* zFRHMhyc0x1POgn=c^8yzs7#}-r{^)c$_Bt?MU{eLT3T8X($Z{&G(2X}pr#}PN)iQ3 z6mi7;O+X2J;^N{!@SFs*>W^ytKwSN?qr01QwQ@K9!69&+fzH4l$)U{(B^Eshf3F|` z0E(~_#)9f17vwjnIs&9xB!~p$*DWVNy6_1L3xhn13>_PQ4Hw1z&J!aIKfw9C06f0IJ zts8bnJvhH+W%+`Z^~aAN%-<4XV#*iZ6;Gg20@NrlL|u=!d#87ll{c?TZ2yP@7L(P> zmlshU+kp*ma$g?mujo{^OsxP76{;2S?i^iRcFh7An4s6F8W`x8uK@;%>{odBGp>7c z2dIz{rpurSFgV&<$Mgj4*^3Xy-Q5-Et{e}YyMhuz*g&Pw2!_%KcWhj(KIH)U9pk@( z-`1tn>$CD1P=zU}YWo0ZzY?Jb%F&kqA@?J7AZkk#nD)x2F9Y|P?DQ10i;puz+gn;L z0)Lg&dfy%RDWrKo$Hqk_Ku4r%LUrj9GpHRnQGv3(92kaOI>k7U+|D2<(>28r=-Mbb zslUIUOgQrW=lF7v`-DbCrRG<Jk(92n_2mjDYqhgvS+BTr$=uwRH8WVZep^M zAwJT*4r>o3CIc^Tahv^}JCy{Svq1uUd}DDABZ{S2s_uufAU~H1p}he4I&6J9l&iv!Jep3PCZj*tC}b0eQK33fE;yN|XoozOvo4qQTzBm?C*v#v>rHhuaRaQd;au~`ROXHTFUH`=KmQeC!7 zl?P%u>0tQu=zD?mnJiFBLite-iYD-*A$WA{)z>UNZEbdwqRTT{hk;3-3uO+dpHY7i zAsp`^yYe}b*AC{)3;r(&FdCrAS$sG?(%Qh5;dB5A0qiy(sG7TWU|;!u`Engd&n-){ zqu~PSdR`~*K`+?>@Fu%@E`sVrp}=_;6B7g*7p2yv9D^PZ&-wi|~|T8IUDY4q4dQCBSW>Uk08n zTV;EfO4JYypA~|16ukHyL>)P>{GpsTsvHoOPe9HDt??xRdvq-SD+XKro|gke;z3B)PK+T^Dzp1ld>U~W|46~QSr zf__5ZzyK3aDxiqzfsf;K{dnhwgHmfOKX9l?e26Y_8Ax0UmeY9HN=08eG(R^t=1nME z1FmrmHpe&!AY^1@YCv@dN_3x4M(L^!>z;<8`;z^1i;qfC4kgsotdRg|1wlxu1=d~l z$pIS>G?Tl9`S>ouOewKHH!>RcI6b5TKtRvPNYz9H_<1)$J`mVy2!G;qVNMzN_MFiD zWdAuTJOCg|nKG=*p<4Hu+`|~|TDU&|e?Bn#%oD(+a~gMx{K7xA0wh2b^#!qL%8q1- zU|U>UY>T$0CM7sAvQ}W)N`uJ*^N`KKGwnlZn8658LSL8K%%X=rz-F$2B3E~B&okx& zuV_r^CO}w#;G$Gx0QZ0geG!(t5nh&vj7)XDyg6FnLvNhGeyX$0Oxrz(XLHyx=8 zbf`QTfo+UM4b}jtMTGG8m>7ZdF{U|A7ND2|GuHq0YkY8bK-jzh%FL7lAr~GL z&_Y$zTOG(vhsMg4q2||EV2l9%CV>{+aWShHW#Ox9X{jB^N1{_Oo%I>ulPdtK0kp~K zaJm3d0i&yJ7~XR>JCr&#Of`KcFT6=)kwkHw{I7&5*Y!(-*dpBL+3HXV-RkGJUB|zXK1*IDV-D_ZW9H7IILcr^`bU~_*e2@BQbI?`onedL8YF8PVsRNR4Fiy~DHuI-BPp(`N*Lk^ zD9S$I(P=&W^^XBUK8CJ7XeZR>9D_2Oov=Pfm7}Wejj3jkWw8R$`wU2|0ma6G?uBPw z9v3j9!_e#7<5bqjccb3H!3nU5(DPz|HUFpqM0KC}h;<`Ci9cEAi_dE91o81%PC+$5 z!e`E$sTmoe1!P-sU;o^GbWe_M2skumx85nECQ^`?0|gn?paT)qD>jyj$6^Q<0a>C` zqFdzE>7i@auh&7cD25A}LVZt9kO9S1meCJiKkkSRY-JTZB>@;SQrH3cVx zLkghVvYInJjPZfv8?Ny1kgu$)811N`Q+3;6I(?3b6~eJ{=TOA9-Hj_}Y;a!8WR+}f z*`g|U^y~sq#HJklGbJC&VVmE>^#Y&ITKF4)Fi4t!<6n!nt@hi7>Dtkm<7>vo#-g~1 z8!MebWNfTQ>9Ug0eZhxwR1^=k>){>8qIBKw0g)Ioq+_+7UO*vHkp<+1AaA39*qIWH zZe~!irz*L^!2{^~3=~_Q))_)h?g_#Jv>OKLHwrIWIXHM2Q6#F^L_wDW)&`DiAO%69 zXJw^Z%7MYkXS-o-n8RYr(d>@haua8CI?j8<5K3P?L3qi6uh3EI7ZJZc>DDr}G` z$D`Q`oWx+!yAzI0PI?2m8}+^*ziWcV5u~T5ena=={|LHQ9vnavIEz@r z*&1-34j21uOP?G@sW#Hpehd!rb-{9bDXihR0X@P(Op@(1)zC-cPSwd-PX)h1u zWo1d=LF6FoKr0JJXXhYrTNwpmyD2LQ>s;x|Jw-@EhH7T`z-YCHF3MI{*U&(Zl7VT* z&vQR^jyP@!Q68q#)zwWY5^~ug?dj=3x$W?+Kxodw$x}WrXHx~0+O;)o&WJ+U9Nv#( zqh?nMYfJ(a_r+hoeu=@gAbiGSVq)Sl>1GZM4_Eo8a8zoK<%cFD94-XCtX6X^s8xKn zY~VlvCg`yEzJ2>PO*J#J1A-*vy%eCJA_7;-UQ9TXy^vFmYT@Cw=*ccnOG`pyr4fw3 z6)ilo%chngu))bFf=(|}f^)viJPcR}Er%mQ(8fY1gb4I1c6N3jI8~&t)aH4qTldt^ z{RlrIF;CVdFz@dWjY7l1lD8(WUv+EK-iF<%(58j{(1pJ>RKO3KN|TDl%2y{Rhn8>- r${2Kc#6SV!i1C>5fA@WkQ$lOLEbhfE`E@vC47o2UFOe;#_v(KE9wO<^