From db4a4e026850e2ad77cb24ea3e0df81e8813636b Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Tue, 13 Aug 2024 12:13:40 -0700 Subject: [PATCH] Task/geospatial tests (#3841) Co-authored-by: dgboss Test suite for geospatial functions introduced in https://github.com/bcgov/wps/pull/3800 --- .../tests/utils/snow_masked_hfi20240810.tif | Bin 0 -> 537771 bytes api/app/tests/utils/test_geospatial.py | 108 ++++++++++++++++++ api/app/utils/geospatial.py | 48 ++++---- 3 files changed, 131 insertions(+), 25 deletions(-) create mode 100644 api/app/tests/utils/snow_masked_hfi20240810.tif create mode 100644 api/app/tests/utils/test_geospatial.py diff --git a/api/app/tests/utils/snow_masked_hfi20240810.tif b/api/app/tests/utils/snow_masked_hfi20240810.tif new file mode 100644 index 0000000000000000000000000000000000000000..7e7189c9351dd87d2d96405a31c285ec1d7b9ec8 GIT binary patch literal 537771 zcmeFaTc~Z@n%_0{Eb=e;ec!ui_{+bn zasSPK>wA~4B>%gTfBpwwTz-q>KPCB(|KNL`@o$skFZ|&1%ik#dACmlgf9~^(%$HxU zFYR~!nwRqL&;6zEUj7Ts@K^sWm&?z7|8n_r|LEoNoBy@T<&Xa1%jK8;iOc1m{~s@x z|Kb0=T>h?K{@LY^{OD(wFaE*LF8|~|_u1vo{(;Xff5#vD?D7Zy>}Qw1^nZVL`RzaS z`Q^X*TR*@2$Zz`m@`IoH{PO32@8_3);6MEQ@}K>epI?6E|N8mmU;WFUU;ekh?u*NB z{0(1R{*!%P4FnIHf1 z^6UPoFE9Vjzy0OqFaF$@mw)6>e|h>=3<^TG-zjOKffBSbXKmV_M z=knM7q3>M&rT_dpm;dFz|IX#_`HSDV{6|0h)#bat{;SJB^IN~V{5QYrtILo6{;w|o z{(t(_w{C|J`-}@VW?q|?{|3ChRe%p`zmf!Ka|FM7kXa0r1={NrLKl@Mo^w0e6 z-|@SU{^_6kJ-_wG{@&m4qrYbf+~Qt=y#mi&0S>O+_6qD3*ekGCV6VVlfxQBI1@;Q; z71%4VS75KeUV*&=dj<9i>=oE6uvcKOz+Qp90(%Aa3hWiwE3j8!ufSe`y#jj$_6qD3 z*ekGCV6VVlfxQBI1@;Q;71%4VSKvji!2W;vevv=Rp447}y#jj$_6qD3*emc&t-uGq z{Cr^8H}w=d?_Pnu0(%Aa3hWiwEAYT8aO1t?hU$SAXg~U+uE59D!DpX+Wc;Hp`+n{h zv;rT;a6QP!B>9-w7qpIhDtiS!dIdi6MEU6PcY?hFFKGq#+3=EnjmMnJr)Tif&)_kC z`Tf*4SKw3P_$kM}xfolySKxzJ;Dc!Xpr8+);tt#^@Vph+Tt07|w%1;PCs*Kg7n1(x z1o&T+cv+4A6#})-)n3-b_H5qr3Vbry^xu&``@HrMWEPpr=Tf7m5S}%;=Xt&5-|&N; zt(N|SVqbaCGamgpF2dJtx_u^Q)gV_Y%kG<0iGFCO$Tz80yKbql2y505)JbMM+513|7loje3NLl>nZCxPVbt3<~6Y-LR zW`Md*&^1Fwh_3%!6uaH2q{LOU3yLACJW9Q$^z;sj-TAgx;MULbZO`G^NQdDMUeuokRfPb*O6rVe}jWReYEtE11%BtK>Fn#mF0d>F=(O?Y%> z2FWn^xY=ipvn!z&X5XK=0tX!$yo8@nbEs(?88hPA+Ku>OqG;4YDUj+f zUtOj0nc<%SyM&2wY6wq0{H;s|eCNd-rA4cwf#z#DqO~#(bQzm*ciSuQO{~CyEHT@q zQe*k`NEcB#$Q0hRoDkw1L;?0x3Y39Vh8Hhy@ub1g@hjz34MZH&%QMR~%Wu(daj(FK zufU-Il*awqa2Fe!u$#7HK`In;WgEA`bh-%G`;U#YfNi|88)N6d4dNU5oGHoeF_MR& zhEK)pGa+m9X|LyB`qR$lYajh4WK{^dh0ErbI_S-tK%hQi%c~Ji|h)< ztp*?yjk`y2&zCp;1I!0~iPF0GH&JwDvfAS8mztZ=k;GKkJl+G zJLI`&=!o03PsqEhGvDY?k)R%dY#mzqU!^%0r86YiJUa@kXq?8yXVap^(urn@W zY;mu^2d+RfHAq*UDGNUFOeoIQemUxy zxEv@ArypfNE(MKwbXGif+bi&qD^Oe&z(56~jTYFWz^gnZnP$1bKCtS7PhNKt;Q(9h zF=L%RwSlJxN|sa%iUe0?PYKszU&qF_!TJT zifrL3GL>^{T`*)A)z&$5mpj8N zkG@W^7`M1r;QdzM34+YDFbk+_S}aroEE-4~$y0^=xkxHz|$+xfI?==4RS}esqL+w z5mRHSDnR@~a8*J_$5_3(g?rdF@CMwtOL8IPbl{K}&DG!{z-P8R#xehCe!II@;O+`E zr+4YL!EDbQ4vKE<`LP;kl|f#j*y5Xc6-?b^Drv&)BaUDeeGhQ+Fjpc8Ry zq&dietsxw^tt8wgstjmoGDgtd^}Ca?4AkA4X34n4Z+!(`%`d?##qHH8^zQ_42hR%+KK*5%_&o7;S1i(Zzp@ zdj;OP0{1IH1Osy+=pvY}Nsa2rUvAB-Vne6n7zBrQYOh4gCB8Wjw4>)N!5U8$LP?IW zw~qNK`{FgArn$fWo4upd!Gu0hEFFY`IZ3>)9mieS-qXlrtgcFu!>G@1tWrHIj-gpgTDm0;gq%nBHlVDe(-6u83Z6O%IE z;O{M)9UBOjXUSIX3_Kt*&(>d$#ErVhMa&9M$wJX*f>) z+AiB#9g3>)dm3YsN`vR2y4>+do=eA$m@T}F6__x*jQMyLf-&>DTMA}(xMk{KuB{i; z#>tTVySC}z46X30FL7+ve*I6C_*C&_77B4`b##`Ersdh<^u_`gcxYs^BZ!Hxzk`#LO2%nhYqcuEVnRE{&_?2|x3G-OEO zI9?(H7l3`#WWmH2yKvowee>K+|k7~;w7OtcgMWBgv zySS4Js?jKv(b>XYfx9cf+xlH223g0e`9fG?@}47~FzyW;X4S-JB`*4-Rh$JmjI@f! zRsLft=2wB%WAmvAJz=_b{yH>0;H11Y zxeGgX^3mEl6n-UfIYDT1&cOWlA?oMfUm+M2Va)O>%So zF`X;nkmbs{6`x`S3ZlAc1}kBoVnP{{7l=8hd-KR#0q@$QiX}*L6@9waabtKDL?xJK zr{Kikjh^#_=BUa|3F%)-kVJ&~z;V6c8)Zha*-36)b|u_el~1dixZRRI(cFN{uO4FU z_7v*x(C}*pPks>FIIE^YDUMdUYmXcUR_D!$V$k4FUI#!*lB?^84CKiq;XYO2DWy2# zaGKHXzM2)N7l%(Y-xKnu?rBnBm~U^(Cx?|q!Omtrn~e&{Uk$O~y|)C8WA_rtWUMd*mv72@O~>Wh@ZTPKf`&IsMK(aO%YJs&64^C z)Jm6vcy~=Ko`IY+#jlH_^LQe&lX-b4m(JC!n2q1cZdkVJ6R$x1z@K;)J$&Ofpl`$R zUVC+WAB~4z_I^|XHU$Bx>JXf>pseVODTC*D0CM@c3Ora?S04+lO2?w}DwOc5KvJJu za-6lh|A|+C-?^tSgu9A1e6su16q$)ITqZPS@T@l61U5OWo{yf2a3Q2b>|(O>UZlQ< z=7?|&ah{~4w8@jb%*n0+`V7H#YBwL^dkUm4|PCa@J;%c^^=qIyrY{+;D=DC zcO}>=V!My%<{Z*$AqcAIEQ{$!5BpAtH#Yfs8VsN`sTZdj7pXI)q1nRAT!H-9xxUO9 z_N>et9Z&W+<6#3-HDJyy72t|+LD2Uf2z7^4`ie=QjjE-t2d7C5C^oKmQU~%$F{Jj! zEMmvSmE*b^x-IM#_=puK!tZBTDRdEJyAG@n%Q?|yAQ|-Gn7y2=wwqh<#fHw!uSdwB z@2Kcq1vuK7f{K2*%9i}Q7?xFevEpHSJ#@asDL*l;a2ixHqa$(Gf%S^xzYbgpBoM82 zvq)DBLmo1fVqFF+QN!iXppX=h;Wa4Q%I;E-{dy_{5r4~;arbqv04o4d!1B~;h`$44 z$3iG|v;D^(ytJ?Y?D0)^nNM^ayY zt{zf$(5$Ks`tP zd(JGQVS-!m24JT#$clh>H*eQvi&q~InDhLMTwA^~RWT_S<>YEI?ead&k&$BZ)DQDd z?tOhxEAX;@q2fx!5=Sg`o4newQMZeo#foHG*fdQ{re8#GE+OS-4#kj!rdy3!(VP)fh zfL6&=hc6hUN+PK+%HWG#;)ObuKpm(OI6IClX98|V{qF4*c-9J-Qj>@K+@Eg*U(?r{ zbE8~A-n=Cn`Z{NXb!x-LSqPR6gUQi{gB&#<&H;%$`4d<(&L+ zz%FS7w07TL-U{@O`srtunH)LNtV+>1K4YFv#Krq`2b5MglhHf7i-Ak{M!VcY|AGdl1@;K69&mkf~UR4Ka~M-f`#tX& zW{NF6_zGlqzP00qYQ59nTp_SMpBQj^rV^9|bXX>tws^4aa&QEc8w$}a$COPJ90`Wr#FLw% zaV82uHdOihx@U23qKtQZq>kA4y#jLuZV-IhK}{%_1+^Oi$(DV5X0}8vG_WGK%7Btk z<3tdqKn|-0o;!Y6(YRw1JG(5lrAks7jB~FFF!4A~*gegC*(;D05Gj*z(Z=p8B71ox zz|xaE;d6Enwe0RQjb-&y2k60NX(+bhs3FllKGcR)Q@x+)6n)pwOa@zmhFTQ+ohPOfgYY_Tn!_I%xh zJ{=@ZCadVNtz(m-&y3|na0oViYdjZHYM@n&yxU%Z_gMkRTASZYR>hU<7mx{k*^hbD zgH~e>qr5!1$-Cvk9d%U(x0Hmr));<$L%B5S#A&)%>FLXquBNrbZE>%_d#!->;V1mF zJ3r--VEB%dcc13pMZuXMD?V=Yd0(pXp>puxL@?r=%2h9bOHC-Pc7P{ky0NoA~g3kv4E5+#to6KXIi%9V|u2r6d~?+T$tx>Z=> zldY1bv#)yvK5zvN+*euLV(aw6ay=TGykBu-GX0bSa6_#LN4Czj!oX&{Q{j`5g?%X% z&x@3sj|PJ5hiiLvhJD{F@S!UJWnt^x8LqgyQji0dQ&4Q~YDK@Rg1&L*2%#R*ssJil zQO-g(8j-G@^w3z^l~BvDZy&t^O_UNX;;A!uKe#Y9fF2wBe*r3yLg4EoLOogB_3ILS z*F#mqNas^z9ZRy}R1j1I1`3C4tJ85FyZ3@t;5GeB+Px1eckH$4ti3LT^7N4)Z2hUT zfU+>r0vfE%va4iN0@1wq)PE`JVRKhb^scmF@r=&X8rCiB6?onX465=n;~@W+n^Rp> zwli_({uq>}d^0)L0K~tDRj?BJFNT~yRSi@?l>wFGXelQ24wagUqLP?a zP1f=b*Oe?A2oIjul|0nOp)N6pAJClQjl81xE{dY#p57RPB~Pge$uR)hUp-weOkHiFJT3)f2@}>guewT?gxGyil2+^9458 z@zlorO7JU;5kAa%xE-{&IvfBdoDmYn`bIo&rbt!5sF38G&G@u=8{lt{WmmcLOA!KP?uByd%uEmVIGn^nI|UVdER`;CxW)=#je0N^i$r4c@{zu z55Hr1%-`$|eS~DaYbGH^V-5Yav+HKmWns6cxgDJuI&uP!M>#BWTgB36LYdsY zzU~#ckJZOM*X_*G+Mr8gbaY%Q3l;qjZaGots+pRRVA2Z^FK;~69X$J0SAH# zfs2%Rzf!rne#T5BZ1SI8`fa_j0=NQeY z&I8GsDOU?rgKC9hNO5Q$>ehIwlpB-Ts=WdqwE{pjh!$^2Hc@+`lx~i@6cg;4W$*X~ zUt^(@+?60nqCSPZ=iA9Zb9PhE`xg`87TifX;LirJ-J|)9*Rb z$e^Sq9>~GIP60zsCA5HnTJ*(^yc@}%k53dtH%%%h%sP<|p6Qxo%X{b_nGb-@1~@hlIEEswr>$BY>a+J~hMPL1C)Lj8Z3g z&r**>bE;HN671_lz5ci#F694YLlz0zJO+VgO?SL)%QSI0n} zl{fp%a*XaJFECdn=_VZ=P*A}0*DulVnt0TdKyyx`pPkdWmVz34g(Qr-J=6+J&Ux&m7tB%%YdIXfo?a_zd@=#-*VWR21`>T!vHt zkaOvy!X-qBVC3TL{=JC*K~VEawT0#A)&4a5;G=$ zG25q-c|6I9q8MVRJ7#^W$wa?G_ieAh2duyio*`XpqUEZaXUSs)4(1}Eipg@=LC>@#6c5CtdXO@)v`T|zq6aRRjt0z(1)w0x^k->8* zpto!MJL);B*w`1tLRW;VqZ}@hLSOuM)-q+lLOvrUQZm}Q9UPuBq3!clSD*npKHqKa z(MLo{(Vi)dc%Y*M*tS&dqlxWv;wmjTV4II;oU2r$aoX zLQ5-yq)vx+I?gBuXF|aXa6>)zw&rS8j za<^Fmr+^wqK?&-}r%1;(rbWk1#S?ow9*(W{+V{Nza|LdnIs?V&*d$s5txyVs?OQye zUZNzz!!+&^I2J}Zv>|cF-YJ)WhyJbw`)c!6ze=L~ zdA)Hvy&It$^%9T66&tY(p&62Rx3{wb?_lftLC*o|pH*wnUE4J9zOG(#<-7Vgu$7mi zZ|$>=FlE?eyG}0Mh@pQD>%NAq>@&%IDN$T`ZR7hM!S~Z{%gb1ShT>VXx%tblZ}#uX zvWzwmhtsJSVjK!}sJKK{GmbZB!N_dWZjx#yzwsp`vbu z-I?p0TT2dBqTETDzd_o*?{9Afo(;$Kqn}tK(^J`G;{}8LU&SzLA;?D07|I%({Ken( zuxRk`HeELtK~;ha&95`32`2gNga>Vx0LdN1Yt@~Z90GS|uYgzJn!LF%*bEa#KmfaD z%LdeJ@%hj@J60J{CN-o|I?9pcr>W+=nalEijbKOA&6P>M_e-3 z8XJD$tY9%d7#!e9)gE42gNZb?sf`fKnZse?S0u&0Yf8?!T*#7HjJY7VvKF^3?iF}y z1ty$lFIkFhSAzMpx0r3F^X7DpY-{O@PxZoVu681yv*+MRWqo)h)XNdcPII`L;#NH3 z3f!mtiPeBdeZd~`X^8IdcY}|Jp{ztw%LnajQeG;fYKUc+>Y7Wx4Qx>4`0tyu7u}L0v{2AH%CSGUH~4gLx$p(5K6Ug8(3GAPYDv& z{yz8hTM6#yihkE8a=qQ%@!sk;x&qJML7)5lBjFJEBUGq^YZb$djT>)$m+H7VpwAOt zVtf=(2nup;yVJAkX7E{AgXKIRv+H1-RcH6z-z)I06*%#mBVYuviTj;}caQpVF316l z{YGMEG^a-9L*3AEs4~b#jm?9=uUvjxI;FYY?P=#uLVs^NZ|z=zJ1cMi^yD*`jsT(Z z;7}0QpiFU$=1fa}7+P>xcc>~?0a^-yWK;vkfj&jVze!lg`(1j1H~U;j;cbQyk+3Znb7#l zptMiMZsS{Dfv@BD?;PCo(j#!gE5*RhT|~Pkyt{)x`+BYK9S*Z`d((zooI@}Lp~lr| zI!~p0@>caoC^@Q;`pK|ue2Xj4Tdkk%(7E%p2cIJ`qB}wM$ssV@L0`bPgNp(Z(Wo;= z8xd0#>PNplAeWp}n7+#6F_&B;6oGSKX>)n)+g^d2D`4JY0ZL=jbh>vN zBx&(VNS8&O1yqC#xH$MaLl@DU5&Y`|*2<<>YC_>-D%15KnfuZ(j;g?wfJ&14W!|^F z0{2$njA&yw^ac8LJe>0Wf@Sb(Z(p)2+o;P(hoaip$cI(Ql9Te`aByLGgbuoW?^Im` z8O~QK;dap0zO5BF{>%dDBU|TlB5&y4nw7OJV*mydd%Bun*q3AldKyt=x{4)oB(gwJ zDK-g6NKugc?05_#c}8Ls`JN5|8&5-BW!Q?>zXD+Pq3u@eK_?Ue%m2n80GbY@8(DSX zHg%vIt+|HJs|MUao0u|(ePycJVp!|0!p{mOLgUVd#N$C+3g?L3e={pEe(2K^-N5eT z@qS&61FqdTBnms91n9kYQwVC>z=1LC-W-E!C+&OHMHPZdfWxH=t;e7DL-G^=nLM&4 zbF73Jwxu_>0tXV_`E;xIDnAtn^tKD~W6XdTUZGWP{%h9;BqB$}+!zFh)lML(lu0a` zW9mRMs$1=#hJW&LcdDAX0^c!_TdBA1b*{hxiLdt&OdJl<7j5+oGuiZst^A84n?L^ z-5RT-_kFLxyH{XB(-M1cBRAvQ>g`yxnqF_nBG=o!eG@3>_U}x(LpZP_O7h>yIRo5v zMe%ZUvlPhZGNfiwE~y3?JYO8rtm)3a?-lr{6=;B+WOjIk@vbhaeS8E;;G86`0|#Dd z$U(ul5a>D$N(zKtod*=+ja(Y2i^p>ySCT54#=F{RbH+opg*UMR-{>zZf(AtFx94Wf zGJ}FyZl+IWZ{$tl+1=>~uPwc6$bAZnf4nYbWaV^ioMsWwanV1f~~%mtrZ1m&f9?D`76pmY%x;#Tj4A+lDTwRzuP_X-?w_=cV(=rfd91*qGj`c=qNupQ=2WI>7YKqL8l z3OI%yf;8Nz3KtC>GDw><0gL2Zs!$%XB+U2V&D-BcZ zTwy-%Ne)vMKJuU%Vn=Llujn$ zq2s=czs0=*&t3r|w0a+L%68j;i_I_zrNI^?$w8RWZC~!@RQlGTqB|I;1f*~}<)Dq! zBaRC303L$~N=HdNHe;x^9Iaft?G=aJSsOBa2svJVJ#i@H6y1)UADOUx@0<8>NkuCs##o<^K0iQV=J+*^SOWqV8wv|Trp_G!ln zwoHB-p;)|0Lv);zDep0doTT&D>56cdgWW35`ml_jU&}H|174+)+$k5$(?Op-6sz*A>>EcH@k--`(+DushxAXg==;4NI8WmaNkl&5k+-xA z0HU3Gb=vZqU4cQ=*Y$)Xgrd<3I}b!9&;U7(CXPK{G1=0ifx=3on`CzNDift1He09@ z4?5~c-nb_NGl2bMlpg3A{?MdEr&RJwAjxkP%^_}g9(n~}=DHJ!Tn1uTdUO!%0s!PXBNQ>*0X*aBy3(rBvL4qVPy}}BOwO}r zJ1aC10^F8CxA=xu;Oh7PuxC@61DK~Qp|syXZjH$XUS5h7zQUD^_O0v%B5;H9Ff$l7 z-Qu2i7LRQ$Fk!W>n72@Z0*Jsmx`1q8#O?=NE*1?spdy^Uc^)+5{XHXl=*ZcQqAM< zO2}gDD=SbO_hKyFMVwv*QoHMZ)&QeIc|#JoFsHBJL<81pU@Nu%)hIw-qhKM(p5&)b z3de+c=fO$gxY|Kc4Cv+=lsr|E+JQSWtnrVoyZ3rm;OqQBN3J0p*QH^y7ChD3b?nDN zj2&8lryx=Y6hR22&o1srB3 z?E}t-KYOq}#`n%1kVP7hGK`_G@g$?AfwmfSClQeQHD?|SDkye+>v1PmJ>pUVy7w%m z)SBa+G{gBI9{qOj0asuFV?2FfIZaa*dP;V?1z zPRE-MJq{R77+$SBGS?B;iIkxSNwQP!F`}(CPJ8bATV8=e?cO6kX`aR@-Ajai;#mGd zpK9=H5N4_EQ<4iLfe8@kCKp20FmaKeI=G}N;!J^*KR%AsMY$6xrCWJ-ce#&27Tp$K z_X^zqwc+S~Gnco_gTh^Sv^4yay=MY(p{D{uzp zxFT>hG@*5LCtPQe95Ui0$iakD4@Ou_JBW=Q^dT2`C-ojiT~~t4;7uY&raU|!WtWLd z2IOjm?jB;9Z;Pd;C>h|cU;FR-1FS&s=?CrJc(IV8bVF6X&?M76*zl_0xj!2_>5Y97 z>RR|$UWOYmULmvt65J;)ftGXI;NTlpX>P(1Y0LhPanz-`PObKDlk;nGx_>E8b3SXID{ zj{P|t%)fgsiVV)#AXWpbywVWxWII31(&0p&Cv@(|G1p9C)ogjMz%@iqG352a;ieNi zy*)(sx_86wI5Gk`WcysmFJS&hQr7}0s0d{w#5LfU9=vLTTPCy~D@l9v%3N%CrvJ?( zTlrR2AS~!<8Y4d#ownp2%wi0imhW7~X3I{#yzU&_`%qA7{i~keek?OC2LROju0+lz z9VOx-^DHS2JIDPn$c@hWC}qyzH+;A1RjoiX``(jY16TQ&=wNY2mxEINVM$sGqBJ?! z!&{Qg+Olb<$D*y=??hwMkFvtJ=}S5Jo0lWXK>itnYt;&_H85H3v{?xy)mV4jxB3mP zz?qU2fwB!3ZM3xllbpesv6bB~E60T*CzgP&&EDb8r6^gQxXk#X1#@R(r>}g#&qvvK zkCDttFdR}7uCxwUq}?k)-8$bW(SPqsnz)J_`c8?ATWp+E?UQOUR6p_MK> zz3&E?w;zUwcjIa)(39OGI0aG3l`HPl3O$I|-AO1ad#FNWqUVy6JlnNq+P3t_D=^@? zUW6MQO((@&g`O3G+%Z%G#o!<%Mqo1^E>s15E3lhtQLcZbt-<5bojw=jPYU>V9grwx zb&Q@*u%*2M%L*8)rGFF;8z_`b*1Z)Uyv<}ZxGla`vU$^Zm+{$vOTh0mGFDo)cM=%q zR8c#1Ii*hq=X7OJYYroxvXGPdPe0^HyoJ32_g7$`A7)1v=DJRzEchUg9e?pO^18q* ztoPPW=ip}YJ6l;pnA=d-6g*=raMn3~-!EMef^PfmkI9$Xr?x0D0Nz-PdmF#zNdciH+H7RJNoUj>qL2>bx8YJ+lh0?P>pmwspwICaQleR4<2&|e)<(KI`5c9IFnN#0W)C-ji~_J-6N5##|kJV zYBuV7ZPX13=qPzIb4{zmLi%j^@^MwrbMDNd1Nz?K&$vG^!DmEom$$J3*A!Tz3&f9C z66oWBgC8}ajx}C}ma=b0g%4OB(ez7=9k?uI(xEtv5b{c)zTbg7BP1!5z3QmEGbbiB z!pAsU*emcsD_{VL6!No2}nbOQ^j) z0sf$%WGaJnm5KXRT!yZ}uII5e54!@P=Jw+YuV!z!2{e~N%0TL}^a0UTFfx|iTedfD zz&S~^t;t^5x8mVfz+By51+a8O zfVW~!lJFbA8pHTW8$lku?XUcG%&Y^m?HGjyYNtpZj}%c?N9j`Z*I?y2CYxlsjuf83 zwS72lVXwfg6&T2f%)EJIr^a`{B_8L?5*DKLo!c*9V^HgRWaCaICsV*UA2M1DHtIr7 zC((sU&_BPd`5noY@^HZ7Ji@58v~ub0zCZX13^3+8+)%Z#=@eLnTr9#jNLR76w&Q$n zItyRvqUpbL)ypX{MgjMF^H)&3jXxqq;~OPG%`tD>OpoMj@WZ_&mqa=WD>wgFBp)JN6z6d| zWaNn8!49@V=iq(ca&X(7x3B{BW1EkDp8Eha{_V+q=m_gj!>9{Y(EE7X@X~m29+f3D zo&_&1(5!CE3=&IRS}o?O23>p{aBy@E66ZtkY^k%^U!)q#11;(;KGq5p2JiFGBP0MC zY~$x7WmPRu2gUACHrbs+_zeODE1Ug_wKC?vIQib5x>E^qBmmi+`3$bc*_Gt+deDPq zclHX@3P9|;7Gq^-$_^{)ke4*|o-#p(#gs*c;{~YwD`QQBL!qk6(`rw)|1=*7!oTW9 z>S-&Npx>RkE3pQx(N4DA-7DY?bH!toub0O-S>0d?FCd*{%PQ5_nAy9fphfA8yhmY2 zuS~9iOlR<#KPLRDMFgJpCCed(_37D!R!+c*WliA zk*(pe8~x$86o}ttWY8|}cYOAA-fgD$2pe@I59}0+j1)ona2xI(&FL(iO@y7j(L|55 z+g^bWSb=lpOxO!pazZJjjolrk)Y2=SKd7q&)=G$vOQqxFMu@&7SL24%#b+gbMx1T* z5PHUa$9p);m@Pf_3QUmA%+19(adV`|tPS7}HgyVMr{yS!J$I$`pA90r!}8yaI5lb_ ziidr!6o4cuNkSZ4P?PWqvbi$&HUr|3SF-Y$L0fv*6$lFhyJxOLZ`kFoB!)DA8oGMr zf#Jq99BVJf%_RWtkZja5lP;WPM!fYnwruK93ze3$1#5+SMytF&kz$hQ;Jtcm`9W78 z$Us%Gc7Gl2(bW;Wtb|To4`8})H)9Bm;n*yvs1#DQLl!5vGoW(#qk1ApU)`O!XM&dM zNq^i;|DB{vcusRB+S22!fYEpk1ZST!Q|UM4>!zzs`{2wb&cPsnAvUfn%Q7DV{A)mT zQ_UMv9Ry0Px2BW7Ucyq^q`{VSjM0g8Ux$30AACRQ?G+Fvd~O>{3p8MO&uJJU6r)5n zf%CvIeej>AuHXJoySm6e7LSL>$ z=|dknh!*bb*rp!i!nWe85k^x0A%klEG{B{-Ltk~!1_dVba7=(pvBZ7bE3mG>Wao;K zK{Eb@)G^e`T*w2pMKZLzHR$d_G=B?#a@B%fl$I$dj54W>onbqVWr~PUe8`Zd)DfKs zJAz%}(2nF_f)%%Kdj*ax&^v0LUZKGfq5J9V&v>OY0& z5zyHgg3_S~igc1HQoYncrs3IG#JQg9^l+zc&u-rzeFX|3ErSs`bJG33mm90*Nc8@u zF3_t1sO+_eN2lwu+sEoZ3zkw)K{;kjia)Mv&qq*CH~EmVr<8m36xsUB3(sz{jLfau z_q_s7tiXWj+DGptu^zEvcMoM1E28xQz4aAmVhWQKjEkm*QI9+et|X{RkdF%?VNT?s z7z8ihEUkx%6*yGemC!5kNcU#r_{>vV(o`HB@AR}EI|8O6_$CxxB#=y_Fd~`c7W-bY z=P|VB^Sq%Z{7}P73$>hX;k=Bf=i)q`chYg5%6ErZntZcz_iEfp@oiUbpoo zvDpe&>Oo&#g?O8Zxu9JqfG0+`Ayhv&v!8>1r7^?V_3ehbN}&3LRVoEpI%Ga{+OvkN z+tKA{+r7820{4I6ljYR)v!lDES1Dwh#j8reY{GAB3X{f&W(x=9Te&>Y>5~6^0$PIJ_$@u^qJU54ZwX9H7v5>Zio=9DB83!yhZ9ARHlbNc0pM z{KMgy&PsK-W#{)RmaCzRswE06e=4iKTy2c(FGG!16*>lEM*dqR-oCn6M#L5#bOjb_ z#tOjAEqyx-u1az6bs-SWsnAWL`7?B&%6Dz2l&g_YT*+1Qicm4Ch`3V_%DkjPPr;Pk z18NHKh!Wk6H>hcJ(7rs{3N!!^R4JJJp;VOQ)Pi4*s!ZTcTB-x1 z{Uy36%-JHUI(@Z^qQ*1cRKQb~%7otvR_v`-L>-5iLEL(xujD~smnk;+)n44~tK3Hs6r zxp#gCl20G9*r)Gnz?L6k1qxpCpa-tKvbgXHv`2uW2jB6qP!7;NmW*u|b)QWYkH3AF{)pNSC#6N+%~~$dRJv>DCa*;OtJ6fH|1m#)OujEk3{s z963-If!afPl#axiu`+bB zjj@?Yr6RBWVq}ulTFBg6i+fw9Zc7ic0^_GQbb9CpMh1TjljORaGZ91_$m3y?8WeYoT1yKbdOZW(+!c<;yPricR5v38Y>Q-fSimii=s z@1UM|s;%OaFjq|Ru40H{TvxB5+`^k&f%9)|e)Oc$m)I=&@o{&^8Yf5V@icBlU3`Gl zIk>4)3CJG3?HD|u67c6OChhsI38FZJ7{{_B)Wx^+wIS2v^@<}SQ()>y_(nDLsPRs` z(@|cu6q?0-dAJn-LPM~yB2lPRQ^%iaZ$P0*&2%D|E^-yH=Cag+jE}ZzL(h&~kj9zA zg})0JS2y=45V}lJX$&nUPgp9CZy$v*DPIXGk14uL+{MT!iuU~CS`qEB?+>y92JG;x z6LZuMJIt^*AYZhNqO3p;V8AVpDFjnsv(C48Res9QF^mf8=94<1Fw80iLQ+N`jV`#+ ztl9U~vG`o9P&%oo7RLv=EkL$iZ>}ku5_m%L zcNuW2Jw0PMiFG6_XWhN?QvjYaUJ>0zPa?jwzNfdZkFWw&0T4SH9C(u*x+58Z23>LG zLIU_kD0t&wnq`+KxKSX=3YM5)PI8f+;l@}N2=zczlyTLC5D!H@ggJejmB+CfD?O`n zNu0;c@DVe#!kFD|54{3mrQqP|Z3@ROT?m+n!AxP@Na-AUc^A zwe674IYAyJV??TkNH^}r)|I01D2f$hhWGS59Rl$rwj4W|6ITVuW7h&}^|%e#4$Z){ z#t2mhkkm}40p$&*XfTAHB2sI)ov;E^`gfmNAgTBYH)q^oK zV7nwzA~gVJeVq;um619!l7UjA@7i>tnawHptrfs-yUlP7GQ!@ zKYk3nkBl9)O?~%PDu_}I!pd-g7i1bRF4PGB$cTrbc6uC5J(c3DlnEVYl4So>fv%Lr zsRAa#Bz8K&Nb7{9u__PZ3|yF&NAlUltJjFw!h@{<8*2#U!EccuIEF)m0@d~xY=B+B z;wlpbfnAnzb5USHS1qQ>CXn%u$&02c->^d0B`1XYQ z{hxd=e7aq+Wq6Dqoaif%%w~WDSZmNC3VHqsEBY|(^ppbB8}I?PBN*D8LTS$mm~P?)AG0jC zY@K3GA0cD`HGLZFdKZsSiGwaqgJ#DFPe&=3IL5Ju)2hd73#gKvRMvYSCgfS6R#epi zWlC4g6`KGDqGE!1NTpNKq_>3xyFJhfK)~|2s|;90Vtfn~z(5e_6^O)3gt)pTG@O}0 z1Ih8g0J?hnb>Wb&8>z#sGGz>>wmohh(goH^XT%L5Q^Z68Fv0v|T<`I&K|dBD3K|j1 zJUFf5O{U%6%nF1dKh7dP@p$nH8la?RU=~=C4O#_uQkm~?4#==v+cysE&GO{xhQgJz zOo?Q2bxxN{vXCZ2)Ii)SjA6=_-9J69juFy`^au6zNJ`<7sEFL`M|PfSsHENP5mumC z*0a`+&kauD=b8spiGlt}m?a2uC755c4~!9~0aW&13 zbn>L|q@9d0lAfo~B=DEU`lQkofpAj#MwBO?m%mn)ZQow`3i$VY_K&(a!YXLg-Ees= z@M*t*iyDXubqIIYcWAT41!9|WdJ8AJzL1I&Sx|TrN1{oMA*s9A$ZubJNUCE4Et6qH zGPL#D$}k``O#v3<Pggp^nn(=%zdaiFz--8;Ke z&U`&Kn@k?K+iifZkkivF*3Lmfmakw?-;!l9( z9Ej{C^Z|gsdiKx4Dki{=b*2)WUkN2v3#jhMoNB?kc(iLn>^fv}0eqyQD)G%V~_nj|7un*<4%5nWQCpO5`fDdetV6G#t7r+WGX^V zPg^m~lpIY9y!aK!VM=2^5HfzI-ZdgcQ--LxIs^zn1V{~C?3rx`A{F4=OT?QQCJxb> z(qXBVXzN>*d$}@Ya;4kZ?IJZVO}o9773klx@AcyzfbUKbK{*&*t%2aEUHpgz(D13W ziiBhT^G@$lP@jqr(w6LthO%(D7g;&{fqYj&F{nt_i>pjJHCGC=QTdiI)!aG5D`*A1m9*OX9HjoWkXH@ ziSCibW=9m?W(-SeL8Pgpc$iCbD#wB~G|B-Jj=IS%QIeH5q9`{{#dgmLM`ojr;z->s z`REa;*^q-C6*{Mb2B+QL#0spxWy^i{10FhPP>2mkj%+uv#8z3}7K+<|Slj99f=cuj zgAQ>#x||BiMwJxvczEz2>8VH;wUHe&P;VHpd!P#UbZee9^`4crO5$guIvAXeD!?{o zLSc}cX|46sx%U07tiTD9@mL4F*$@W{3Ko~XSP?9Q3mOW=r6O_om~i%7>?RiR%ceaitf~$wI@+8+jbnWgcbBjVo&V`(CI=QKK@m%%E~GSPj4sT7-C3ntgjy zD{$w>4Fmn8=Ak>f1BGq@Uyw*VGg813q#( zbt8~e-o49&8+zxkE+ZhS;csDF=;rA&ccN*`Pm)bckJ-3McS>RZeR;hr!0$PJQ$uy* zS*8PgO@H$sHVk5It$3_&=AK?j zeRG_K$FE(|wL6KNog1TjXcvaP_vIC@Kz{H%z$C_fadqp#z)#IrP(S&&Z)Ua4;}^;+ zM0Mf=!AcOSWQDoCI%CR!ktiAGs+1_GDGzM{2gt9@I(bYg6TneM-9Xp4*Tos`T^wfe z+k^0`0geVY34Yx6+4tAF0{q@r?Nk(P4^8w09b5rf7(;1V(NaWoj0j$8^xY&72Ekw% zM@xDWLcy1`1o5q0SmPio?^l(DJeACql5nS?g!)mKi{zwqRn$O?#eg&>=}zVYRs7`Rv5a% zK4e2JuZK(>=N(gX(2|2<`h#5ISO$K)Kc;Xl?^@QV|*0C=y5rT6xex$b=ZjW(JWH zGpjL^!sGOXxyCP+tPI#|0oAswrPtYtc6MBPI= zVRx_%S*hu^vt#N}EVBTa=m^=ut6YJ`0fc4?e6v5KGcq!CBySzcrW=wsXTlnjLf&@J zB@&1nuDeMmKsj6`x};KYrQsyE>S`u($dKfrtxHW7f+xkO%OEFPpxxu)C5es)AdwSWHFMZ`UtZ@5=HuqP-u-mN$HT&cVFd~8F{A2%%B>Ph(lqC@vYUJoJ39!EW5qN707RR`CoMpOnrP(h9ocD_L0#*q&vCo|x>kT+_WIX<>KRt=Q*_XP$d|*hXLTMV|eBnxQGhzP+XuD7-#=%7>rd z)K#tpU0A{|r-S&@$6nGq!!8Gl!T&jD(sc=bh{qPmZREvqpCBj&yVIV^7$Fbm}0ljg7HjV_7YcK;^n6?@;-EdW&TmvC`Wp;YM_>a8l zon}WgmEfR^n}qT6R6JYFju~V^j5rdNzZFy4@xyHuZu+hcTYhaTu>A7JUGL6=7&rxi zRB;beAtM6&q!qaM3&Lg)yj3pJg08Mej@5A3{3U)&MrLE~Fa@`_kC{11sYyo46pTdj z(UQBTgXt15 zkNn+`s1S-7vdKcpkkdF)9QGWmS~!z&G4#}x!FfBI?rma@t~sWYVI>VKW0xG5_??L_ ze%aT@?_bo6z7xpZ5p@cLLimguF=~J*sZRw2t4(Yta6T*&TnWpmU~i?cIWF>>+qpt2 zt_~ZI=)p6FQO1h`&Wt^tqb%*qnz*TJp>lYK?@;HF6&sP<915zO3Up?m zZ~d7uQsD1_usSvwbR7vgQQdp=le&YSh`fr}idVG)V&py+hbJ*iD`b-n;{&;)5E`f; zSb0QWAzW$^4XH8fhGbobuM@~j@zzWzs|qW~QOG!|tIQt3MOf8B)j5artOH!h<#y2D ztDu9^?!3|!XsDhw{ka|E?$E@8B2*SD?nxs>(<2Dkw^8Qi(?qOxCPS?RoEPY3b6;jq z_YkXte3i8|HE_o|JEw&e8NulfBShaZd93ao99wz5_vmN+?!S%|c-D{Hk8(uJp+^+} zHTGES($Tv29|Ll-37>Mq7o8~gXGY0_A=|pLi(;3$O_6y17IS>TFanUKDmj?dfRusP#9wM zmnHF?(s2oKo{(q85-BHp)`%l70>wjqX3;aL8J3@?eN>gDZyneDeSd>1P(SiBo@{~d zhPO!?+cCZ{3!%;b!nBaQ%4S~SmQeE>0nSZ78C@7j>p5ih6au{KBF`h!*qn%@m3axZ z=*!iFa6Y$d&Dmr3U-b$Yl505cJrjsuHLH@S5fYPN9OnWxaoOtIf=*5#%uWA#1WGCN z`5+i&WxA7>U*nq+r0U9?ZBk`I*gUF zKg?RB3UHt=+&LrcFw%Q*?Tsf5+MSe}D?xY9rgD<`Fr;(jvd#oQhb(rBC#SP{_g~=( zjNfz!&wGfVxA;0<=nYIkSW121zV`TVY#(Cnfo>otQ=**DhO&DsUYk6KEIozFE`6U1 z(jbR_+=*S&B)KQ&wC79P&-&fpE8rCnATV-=69`uG4xg72Lk;e*=&c>mG@D~}XpV9@ zL`y$f;jB1yyrXAe=gAHa?Bpdwml~4K0iQfL_-G=7(~!85A0NbSzA*l7kG%qh@|{a? z@EYut*2qC2?dj(j5iN|3*tUhRt9ozdr0BeY*51EFLyE;}HTgcaNDE(z< zywgZwH!H{>d((=JC{z@H`qY91Lgb|n0ra^zdbxQNUurZ*4Eai6yxksQ1wP=HUPFK` z%fzj=R~Q8LP%g=+dhZUS@Ery-pWZwZ=$@7C?hW3E7mI6ZriJSlH?N4y0G>LT|158?wo%!yNw6O*v0d{7C) zSL4ZG(ji#BYK3kl0r*anY@85i)SOm{T@pkyr157gCgypiqpXDJmbAXwZRv0EkyccyrB<1tv0g-XpijK_~pkrCYmKt1I`1&z|RP3B(5jH zp2i&bs4*dC6oZrba(wz@FHgaBwwlINZ0U`zz+wA6Kcj(Kl?xY!4jPL!;FQ@_;7!{I z3tk!IVs0%6X?5Xx$xLE{yXY2C(N3YeBz={Gp-G|iHHW0oRgG&G0!73C%d^(lXw~YL z&f&Vw@nG1Uhgktv!3{)O1%9ELxFH3Ay5`sT4 zQxepI>%arUIa6?nDUDLm<3lPg`no#kzjt$%HR?*S8n^svSKyt$t9RZ!-qjZABS4^e zhlfXF+o+~TAq+eXw0cFT;#?sviQ6+t#-yMKP19>Cvq zKNMVdZ8mVBIXu_z645vKKUwqdbWMr4O&EX7EnQ%|XdFYbZX&g5aSF)#X#J{^e z&+FZJw#T==0`sSvQcdDYahO~Y+F-3V1?tld;a_};OTea5q{G2i;VJl25RyVQ=}2QW zBuh`>peXd6kPJ>!?oMU7CPt|F&wG(_s!8^Jxm+ox&0{!f>MT7GbMoBXSG)r3$A@MX z51O?zosYYA1I*o+8g|1=1ps#VF6Kizo@96RM9g5+fys6)Ws1vyE)?feY*h^*ISojd zw3+P^nCyLB*>9Q2?DjA#Fn@L)nVL+ckjo@yKdX?;#c4R#x){i z{?9ZO7`>6R@e5$U_NAyND_9$E+4$$C%7{>vqE2c+cPD7a_PTQ?3k9HSHzFmYw7w$` z55{EM(03^opMJ8o?)zI=f#pYj`jPGBuFq+5a_En71M^DG+2N zBK6HD4&z#}qrZiRUV&`Sg<7yWuE&+mMZ(JwuRyRU&Lq_4ZpAS}^$s+StiqMS%`C?Y z^Xc-aSz#Sb__`+>pOg*D14U>1Dybw-z^ZR#?PMGGU>J zXo7Bm+H_w50A5e62z6j++shdaG2YYdqf0hyV!+#7!<&sCJuq8OCXm#)2bonCRYnjAiUD7PtExz&<7k1IN9HwW(J~NI83D;+i z!lh3l*I|(^b*ezuZcL3vBK}bqHQkt;i^HAgD*_i+g+VFD0w3bQ*Vj8q6phEI(Zo2ET>j<-NItF(V9z&G4*Sb2Br0c& zNpu|h-A{J(uigiS-ORO>FLDL0cV(elzx(;;zeUz~D8CqSCIDiz-o1TB46oCsOFgAp zrib{aiE(by%Yg!;k2m!p(s9T*x{{{`?fa7Oqs2=7(Rp@)PrCw7|BhIX^gIc?A7x|| zkWmOk<67`qiOI(USXYyfN-_Vfjg!NTQcIKC@lvFC>B?l~@v`cu zrCE>dzIZI+G0X0~u@yLguAkj#9G~&v(xx*Z-kHN5n8NMBH9Z;D>M(@y9NrzJf&qPj z%|{nAFKeN%n>6u8MI~^{-E~1xlt*WM$Ng#+ihsWh?-*d4U)&1t3xCY7e4h{Si&Zz? zVt95#i#C8y;^biWsZqE(`cMg!4w9*-CuiM^d?G!Z73J=Hp279=WnS~gIQORoiMj3e z2rKZM-x;5qen$6Xfzn<$_I~#CoDEYC-_f-vW%7Q*q@tZv*rWnXR>TQ%)dYTLfQ-W53ink&KMsg6#5cNw--XhyB; z4jXiF%0_LH3v>Jew$;4ZkFGXynyzDnJ~vw7eVz^@@GtYM;e0S?uDl!J;ylME=|tl` z#~P|F>=k%w1q!;8%K;>H(x3DNZg2TPepG}j!N^KdR^DgaDn|xVDO9Pf1H051q-a<9 z7=M?7B-g>P%I`4EHPe>&3d|J~o{b+gQIpCHaV+Bm=;Zbj}aj?5n zl_1`(8~CDyDPiw{U+6uvf6A2TPi%!+ZK8F$~leIZoIX}OZL6SmWH?u?RY zJjZ+cyulTi?A>~}z;P217`|%uG?Uq$DTz>OS45gl5cEK7p$bLv$T%gyR$h4m{BWp* z`8&oI|BZlzmp)3dWH|($(qQu-ak>{3t7C6He9Zg$dRE{XgP}PwF!ZN~Ykl8iuAU&7 zL~|Szo7vW>ftU>wa`lyDN25Im8EAH)AcbUMSPb6{9!9l1dyOtO({Lw8_d02Aq zRL_;XCWVp0sq2_(dF|U{tw31Ko5gtYGeCJ~F)C-%|KHwqXi2W4&~x_xKYF=j6dbB6 z&ZBn5LU0K}V9V7#Z(Roj-@Ah&YM9&62j$7~(N>f;0wO0kHlYj58YuLS^#wmi*fhF> zHP9${@^GiwV*}Pm)$N$Q`y+?Um;1V{@dFv~y>)_uPBh3gwH2*PD+5uwK;w3e`n4mA z^d8$uMpDtM3n<{jQ|%xrcH$&QL2^$B9Nhs>2-T_%JHDaLqzkn#4i0SF$By7(m*|ISs4YrCGUI9EcZx!|a# zs8qP#aX)wKxT}})ubK(RmCc&0FGa5PCmq2pmTYUKhci^Ze$ZGl3m!L*0ccE;#NCj$ zz#ALQ{%8hr>6Cluw71Lve3D=RRKoa!N?zR`Pkk@-{$Q2Zzfv%domE)sgyI6pd#`*l9CqzUH~x1?bvnX-pnu*DHS_lzFIDD z?E`04f+2cg2iyW=BCL=C8Uu{NoTe4PhC;&^B`me1+`(=QQD;{$OT}h>X9i+Jqc+C9 z4veoM)%WU=z_rZaLCP=9Ze=9MDr#IlI7ZP27pm4N22?in3i?foJd6Q5;Z4ksJDbEP zHcZRLkQ;?84Qnz{k3XFOtN*c@sFZxe4U-n?GVtg&UrP^90x2duC0ISxmdLo=lSL&N z$1+#|fh!u1DO0fKP&ZQetXqML(8GjLx5sJ0w=J~7w`M?#sD`&MBPu|vL8uI0J_>8J zX4-E!zJ$z12)ec!QFt?NUp&PM@A}O@gyOvfkN#`)XNM%DZVyyAp>a>gJeo> z{>n_P0S--t!*}lNhU(a+xPJne&QmPXR)~Qt@~xjs6{Sxf+~Cs;ns{HuFEEj~2>ANIMY*cU}o${t@Q z(Z_gtEUWR>SW*upT(lTffxs~IaP<+M|2^Zx~GP~(E~Tc&W~?~qV$ z@)<4NoR>IPy0p-`))i>rg+fpRq^G};0c;>|;#jh%LzHz|hSG8;>3)QKxy#)B~P7>vVvU7u&DOW0R{sS#8Zt4azD5LQaiF6LPL9tRMFd(X# zotBkILqxE~8OMcEf`mWM18yA@Zvr*>gBifN-ov(vL?z}A^Qr21^W{**Vnzj&sO_zi zs)bttj|8@WZXhGLQgx~x4dMiyps)xpLJFVGW7vZ4fEBUAADyNu-XjCyjZrZxug&|? z4ETnk0c0O-l!#HnBH%6eWuR_xQ$`_!c~nDqH$WHAKx*<_APRrLknTW3FfAk7OwqaH z1#uC!&^%!cg{_GYB~6z3w(-Pbi@z`fx|O!uC{`3U&O(cAB^bD>M1(SIVJr}-v4#;E zROMtH1A zx$wC8>?Ndth;ix9uV(<8n71-EH7H(VIrm7(`q~u=JBK${Q^13VDjpxcvEnKd3Wj<^ zBC=a9I)VP0vW?DA@0brU8#FqlwphpyR41e6n4KX=W4O7)P#o37mt}yqO=?bdG@Ue* zjd*#p`Tdn=C<31;CH0Whfn0v93%Fpc$hi{;F2^$vp^IZwVT6by%kgKu1XD5@t-`d{ zDwt-)1Wl-KW$6R*y2U}Ur=Q4xZ7Y{GomF&Gv9)!3WoB@y%-Co?%W&{q7obGwXcoqW z?)uyefXTyv74axw!x&PdvBA(5%AOH5N|NoI9!bWof(kIvHGE1 zwJwNzK8HhWem$VN=lbOm)Qf2o@VEyUfclry5m}TEF^J$|Ws`)%$h|@1{$mXh2-6SF zODQZ}moG~%-cn);zb9Xq0lBAH5pp%+*+7xbRYujyOm33cWv1{n0Wq@2CBp_+J?{EC zVI3)9B1I~CN4+4PO68-;Mu#>pKw^$_!`&Fo7lOW_*3d7_fNvt#kzBB|&bPR6Wa0W@ ziB=NVvmr2Ag>g5DAY0-(&nIi(-J2lMX$LheM$0vOD1yZ@)JKI$)CXoNJk6kUFW98q zhVe{9pz%#1*#Mb54HBIW zCfA_gp4toOMA0&CdYeo%R-V{8I*bgyD(^ukMht7=mmSs z0L~R90y>8+n~hC`6`Mv3#))x~O@%ulLirMcc&Vdt<={9|lc5 z-99!tfo$rlGZ009ZMLnp(L*l`x@C(XZ9Y|Ynz0z=&KCKR0Gc`ChQRbk%&KuRkv=BG zr~r>fh&6WrU72Dz^Vw>lnN~>7Ij&uBr2Vw`7LhdT1sT8#tL&;(gWR@vuCT#T61v^N zsA^pM5c26K>EX`!Lz{4dk)dQdvOy@ZfQKOXRI&xJ7^7{hr`lcwYKx`>W4EnkZgVGMf3(Bn0^K}`JEyOjx^r1@Rz4z^PLc$;Jx3MPD+7yl4 zlTq(n=2W}1YH?gE9AoIhM-Mzy&_kdmq2mn6`phfKvmZvBlWDFLE>E0-{XJCU`DReWevhp)L z+=#J-^q2NrH`DIje2?~}#7=YF6W=PsM&bkr1#Y@1bw>+Ro&|7BBL0KM;SY{Trnwyv zxP207EFO3K0$NY5t1^-<{Ng_>a&5c0-ge!dIO_T@JP@E9X5i26nke-cNom`9d*W5Ma<;n2SWcct}j2A0lN`%;JVK|q^gB% ztiW0lwrs6^LQ6N>i5--@UUWk=^;0;qz=NFM;36mvsLIoqNNajeN zQ#0l#G@AN}3`ikx)l;qVGXbp4zNYTDSwVc;=NE!iTG%z-nW+`{x&p0v$8n^9p{|F1 z9g;z;*DY-pukifDlE)rD%}D83FL`V7OEVyaLAN;R1d0+Rx!TyFu*UW5OF@oD^d?~E zE_8-JAfB4WveQfvbLA(@Ur3774|UNB`!xK;SL?ny16JK~u#wt9 zITpV*vr)v?Dh$PsPC%>4ccj)1=7K686@`Xi!|ufZXJa|yWw2VZ+8hiUdt9wrC}^L6 zNG%<%qZM)lP5(>=XlME6Y+J2w&8XVt8Vpy3Yp!G16kkyBs2P7C5Qe`)<#Q3}OtTUPWq*}T*8>1TM*Q^|&I6~Yl;-SFsvaTZoE&I9*Xknui zDK@Ja+fkcG*clbvOB8;1Kwc)B)F{Nm8Ht0ay*5}Ll!xg6zC9hx@l<)`Oov>d8PeO1 zw)kM-b_A{P$qewmY5=$=r&_g+_(OPH8#n<`!K2vCq}}+@9r(~3fzXF~qSW0z=VrfZ ziRHDdBFEXFqcIYS0#1SAMM%3rECg`0*z>Q<0F{XBHL8?8az61;hdOY6>XzUgQDg`Q zAUKgH0CxhPo2(Vk0%U?=G#pvW!fl8!kby@AZ^-OGMB^VC_J+fTzAXb%f8vfxby+>A zVY!Vj#rs&TaKtQsMT7?r?x_41ni7CE3J*7&Y@mk?qTy~-H6jH~gY+vv@q}YbxNOol z8EM8hW*`azw`_c90C}aFL}7R3Q;!JIM>ru-Gw>BNQ{#sPKQQ3fqA=+^%2~wuu>fOH zpBc&Hj40wdxeCVOYz^{wAg5QAtIM8U$-s+qTkB4oxl1PyXqlWxPH?(IN<~t^X!X*4 zlu8}kC?mr0>jEVFIX6BaZiY5!+jMPZ&4rFHPbOxD$J&WS)RQmGK-p<7(b^n0!_MYV z*4heu7?SeG+c7@4p$Y{QwSnDWZ+d(_Jg;G4-GWA*X6VdWm}u75W*}9+Z1lxQZQpGi zYQt>tnFuDi`k9-b3GyhQ>qW1BfIZIz9rMoOI0c*;X!*M{Aa8nBC71=Ru^3|xE$w%` zv5N|B^>{%>oy%<>4EgXw0VgIq8O_%J&OsOvW2_^1vyVW{GU3pIBuT|gne0@X|HT>L zou^uzTcFlNooP+@x=1@pflFVCUn^Z3z~M%aQF|&FYd=!opa`G&iJ0$LMYdG)x#8Ri zqOs4(K;A*xRh}vjZ?XKja%52*(y-d#+q`~FgjWEn8vE1(6uh9O=#M~rWjQ!!UM4j9 zOL&f9Kwk46$&rPV)-}b!TTefkfx2Z}i5d=eRZE&d6}jAOgFq8jaQKn0OQpexDV=~P zg|!P|wBMqlhKvAm>tlh8N^3Ja^r$_~VU2XxMg#n)`LNK`TQZP0N@%V~Dp}s=EY^4i zmQ8J>YM0W*XpMl9(fhr00o(Tla#5K905}E~pB0+hI472-(EY0qb0drugn0hc%<^PV ztYl%+Uz>rnZxUl-7-?f&gm6P9l++ttZEEF&s=?0;_`%@yXXQy!)IVOEr4jfvRUKZN zML;+477Emqw%B);5G-V8KC|hat>%9;19=N&SGiavCe(qQV!#2|6nJkT$6tomGYpq> z!ZU=BD0kGn)<%>FH-a_+B0-VC_!B?jY-w~9A2?SLtzuMyWrQvs4ooaVXy#XBK=&{0 zp4BVDBQ|zDNeexQpt1Q#D-0DRW>G@+Vh!!=Y-8<->>AREZx)|m&SU$tfpZs)cr*jF zIkG>?U0N`MRZD3dfZY1cN5-(W9K7jaz|VuzXTlYa*~b++)^i2D#$2;;4Yl^r+*f9x z?hOCbNNdp~Hl9bZL|WmJu*&^pg|}=H;%xzkj1f2ZHG9IwkM0J*I5jzcnN>p=p(m~# zeq`lH1_Nhuchbb?GGH~s5uDzeHO|=CpjO+BR{*HM3H{{)W``^w89xzFDe{bn(|8Fh z*y{g$sQQFYfu3t;gj$*mX3QD=u=6v+|GxZBdB2#Rq3qe^wJnn$8>kysS;VY`F0GKm z;zJQV7pxUV#KKlqzeP)cWYHrX)vaNTUxB7sxrV(q{zdnihkq9JuGGN~{-X4<&!B@^ zCOZx)E1*;nU8kR#p%?HdkbWfHq>fo`$T z1TI&ll_rhgv1{`$EJGH3Q!C6TOp4ug=%W`>aV^Q%i0JvxXMlECtAXx&b6t1=94$0# zj41xWl>Kx{Z0V1dqZLTUZn>_0$*rRY;!t@IV$G9x*ej_Rl-OSZ3=#D7`!nGCm$u(| zi|~(R4o6ilfchpuW#pAaC`{?sl|;D?%!47vyNiif_?s|S7+379ZV?OH%cyOo3w*Oa znE}~&T&)CbWX`ae!*mEdG}kzNAk13SHvTP?tA;gFaPFXdG(FdG7Xe*$@Z@A)&2Dkh zX~xfFz^dVJGgWwqp%t60Ve$U;k`n-MAE56xT4v-#yo_nqagO6b_@t2L*#L?n6K-XP zPpVD+NCr|7c*M5RySG7;yFw*YEm_v@PLrOUFz6hGpD$m!7v}=39K8^gaRFN)<{MnZB16RtGha^D4UC8lWTVQDQx`^>`bY-m9W{;=Wz?DOvWeHx zy3`5Ar9+Ya2C8?$^z=#f6VBJ(gj-z=$dwQ-lNV?hZ1!jpwp%9ZX|SZ+I;`|6qG#Wd zf$4@d-JNweHyeL(ujxF*4+Q;+312l0Z=oCef6mhwV4YJ3%uPT(r?iZ5NVwJ5h_U%K zl%ackX9jd{rY)Kn7SyI(4$IezWwy${hzi@VduYT0_{&)R>&(5%ykg3~Ub^XeY44_u zhQBKV!@lcLdE};Ob|S4O6Wi#T!;PQ$`km2)b)}~f?!pN%ou^7?7J*iA*RbrHFYhQ5 z50;GG<9B9Y+CK}E&M>x>HbRtXgC_pc9O>W;G|9xl&SJKHLLZw@zdI1~>9UefAyXg& zZUs3Bx5ZzX0oh4xuWj*uZ*lu06#p^wrA0x5RSnJW2V%k-&SW}Hzq&feeR>i}UQVW= zKhE8z64^9k`RlHPN(K^{08Rd824d%_K4W5wVGb^{$RL}M6yy)fnM4QRNf8wk)W*t# zY;rN1re1`9UJ6m`54EprIo+nW1(B~)+}oO&kj}uDGHO%1D;wF@*IweL zmB{u+8O_}gh}3w1N@Upp(X4>c-{&34GB*ed{Fyl>>G2mckj^huv$VT7-KT{Rou^5% zC)tEm8jT=v2UV6y_6D9VvA%@im3v*H1nL0(H`la&XJ}m7+)t0P4SpyCyuotOy{0|4 za*deAIta;{%6P3nHY%2JwMx3(Vrz#RZdgwdnQPJsc{&c2MoHj}k7oZ18IVfQdf;LK zovP;gB8yU*7%Als+ro(=Og`jf|Lp36L;Ou&i-k{)gD%pf;Oiq&!nANh?5YNEMiW1u z0ohaA?CB(6)!^%Vs0Qfp&EG+(^0E`yLw=cdmk$vZr;ThxCWQrl`nH#4nhn8mu#-xT z=}>qGysVk6H2s1M3_F2e4;reJ8S&;YynISX)rvKhhRYgsrzbZ4N9LC{7!3&>Rq%P% zVktHYFJ6g+qkLK%C~`>y-wCR7g_VGq7k7;nq~Lk7Kt*C_cynH!fwCjE7u!>b-p(be z6Ewt8A5w!{R%inVrGjFGC6r*VDFfTMM$Oy-&(q)`WR96w{3k>l?lv*eVzbL7XShW$ zCOKQ`H1kavu#M2|wrMJ7XCWfgG8n_2aCvXy@W%?XeH<6f6nXPPpKP1D8Bgri}dpF|*w zA9Y>H0vQ24e{Tj#ZQ{-tU7#gH@yPig=0RDo%n5c-17H>?Al#b9Q%A46WL0<+5lWR7 zP88Y#Nl0Px+YNDm&Q=m6*}O-RJz94^9qwiA`Dko?-Os$!S;qHGu9%Z3PH+03yaCp= z(u|j@Ys;C)K%pfu@oY)=i-U}myUdJ39}pKF#QqFv&U-UpJ5T#zBnYK>u&r|W6yg4+ z9t8oTdkd4v%M6R|q5_BN%Yrf+kQML?%`V6t|Jd;|JZ9R{XNfpA+S1o&0Nc%Hf^c=A z0`{F?EoJ5|7(;EoqB_K;bfc7_9tZVJh`#Bi=r{oy?eCzdFVCSNX|o*GGKco-lD96+Qd;69fUUW*9XFEV2&UCKKN`vxJy@Hx*$*J zF^E{T?G=(L)s==&#okhCMrIA?Y;PcrXyTO_*pxL@!D>gR-CC*7JpQ_cR*?@FFu!s! zEtbqJp9T|f!H`hqfZ$u4=Ma3+$|Sd>ilC`AWDP8#i~g4OufmN0mu@jD<2ed6I{e&;i zf~9OV=vB_R(Vm~paN5*MGT?jpa5bESK0{wvALSM5b9dec+bO|LkJu0S=9iMkLVhjnoO; z%9p{(mmhK>k~1nj(Olb>bUQO=&RYIR2I|SESV+59_`}3D)Y!h$#-Q<`&1R+1dNR6c zkAj#c+x3z+J{aXt#kuyt6bUpY`-cpH+p}jfc&OXL@5sPr%W)&MjkQ2;Z1BdaS@KHT zncVTi08s_BY6R(cT8Z><35q6IAJhupJS>?eEF+MHmXJzd1>2KPXJ9Pqvv|^p!J(^K zb0HB{4VhPH{dW>ab#iK>7_5_7*Ek$GmQX0_&gHUqq& zd^wKe3S@)XuvZzKodKB{K%@@p%~- z6My%KKvkl}P%M>6j+O*lS;hetRVI>gUh(&5HIQsQOEe>3;kT2{;?HI{B({|e3U}h@ zaxI;V>dOZ+u->?XfR@LGCG*Hgk|6@Fhf-90pPtJn={0!Oek%{H5UZw7oDguvBV7bD4KaG)ZUyft2|L@);ipDeP-Xsk`Ak)gLskazr> zanmN}_U7jPoWMmahvE}xTH$r1DxFj**}G&?4krvXyglFG;lztL zS%chMTII!{ibhX9k^$Skz4IL+oh0ZU%P1?Ujj?qB5m;R;J#CdXPHwW*e2sh=2n=ro zVk8TR&QlI*V%!Y1;eX3O&HoeY$~K#RwG*Ft=cl-HC`xD6gulhOXo>qbVKpK88Y_a@ z7Go4ic@)_k2KGjtxCiF2$6Fb=vW0gH>d&ZDh4!~NC=)ZAiXf3D^>L}|(n2`Yi6YiX z2$_nh%JExs2w%R=N>4kMmrzHcyAz|ybRr$p#1Cg+yY&|i)d$Po*2zM!U>#UK5lPJh z6qIy`rCOL}1mRLsJp@aS5EWirU1 zdk{gG-uLBT-hSUgRi)8&83J8@&6-_ZrAOCS^&0o04ETmPv3HaqnOVws(lGPKo8t>VQGeDX5 zCHD=?vh9RsPZGJpzF(P!3{@#%>j(w-nn?=wsSdQ%we(T#+{;$>$rruJS`aiY012B(eI;KgusO2{m2@1- z7MlJ<2BzIP2*SR0sgS(2RdQ|iSN&^V1cma<{xFx4<*SK!q|4MKf*n1sxq{lS{96@vtJi5h}ZAi^yM~BT>m}Q8=Ew>%e&&MN79Cm`i~;zWrOZzX+nH z&aBMAtu*L%ZB<*um#idKt{6734)Ex=QH&zP7SoU$V6awd4&Mh_6V@8SRbi=VD*SeS`MM z6FaFQ)&O(&ym9g~w0A0Jg@y(}F$spzX&uW*(|1fodkJJ@&W530f^XLu@L<@8tMwtCdjra8#6!VJiCRuBgR$EHask&q5kQ|R&48CWD<_M8L@ znS?X{9=}0t)Uead@7_X6-Q82oX!W11rI0pXc;dg!VJ$+&l0r6Z*~zwxe9YgFngQF> z@5zAdIPceH(|ptEEqVhf%ee$98a5fNk8{E;Sf+h8dTq;sR+faPS2`TcFBBR;Jve!W zB!n_uxmiYf8MxGL=A#+d!sbIU)`(j5`4am;crnEr5k#&1Pda8o=b{8{_aySwS+q>95Z~Y++EFpvMNUOJ9+Jb@tcXTnuw8JDnMM!g-FN z!HRx=%3Q8tEp%dNNs9umLUELdEPN7*S^sL z)q4tqTzE(a?j{T$X27Mt7-%{8FZeNzFJX3CsrFJJ+_7l-)DfHiwATm^#iqU?17&}7 zVn38k@l6FXF%tM5St2N^dr5(vOa?1FMf7zA14?Y>W0PqdoP5{@I%jEY4)Y>xR&36D zG7uY7i$%SM`1#?e?c$?nJ?65Td_QqHzAT_rmURYOq`5{Z|HWbRC!piO>018?1~D5O}p20-`dQ98~; z=N2`r5L5Z<0-nYJ%lWe#12N=&d zkFzGp034>iJ`grxI^K&_#ia56;ET>Ju?ry7hY&D56fW@1`lbx{O_B>|R{%TyA6JsT zz$6cs3u(UV=Qt!*F{W1;>}O$fc-*xF16I z(gdIXk3wlo7~J(S<(oe9lMpHeO^o?z1txzuO=4`)GE9$Oo`JgkpAFVxv6`BqF~B4M zWBS^&qDh>K+a7xZh%FN@J(zO)EFB9|4OM$b_{+!xBQ1Jy2JY_aGs5BfsvD{^XWu7i zDPHvdi=)t7hS5l~Oxy~8l$!PXp|!@BW*~2=zbEElw`De40-ikNadGJ0nSn26sGH|j@;}(tcM_q5+8hcZpYj5<-7{cPmK(=9a z-N3{kN#zYDU2`U?-YjpF{Xzyd8=ub>k1ivb23hX4NwzW1@PGxi!5bUz3FeAd!J7Ehie!QSfl*Bo zCq$hTfq6>fk=z_5c$06*fbNQI4o0J|`d{H!JTIrEVb{pTIU@x|SI5tokyERy%-{?k z3?{&fL+J!Y)XU&~m{zRKE`bc(<7YDP#7%mFsg()l6+xNYOHGkHZH&;FF9W7Ge-bl8SrdT;f!GB`6g ztcvUi5m>Tx6_}MkpXvqRh7YYe)x|8#qzESF$gt|QG3nkTN5l8D8@I~W2@QE&cO1Hw%eW^e-YFT z=t6tXs(k_6vpo4fw%I#Zcaz7tiAIRb$5t3o$DUjIvJCunBmSy?e_?~iVS3-*s4)IwgIO_?ne-fp-_|B^ zmMMB*Jt4FQR?{j^&A_myOWRj4OT`D7?tx{F_!Xx4i(u->vYmm*vI zHt_4MVFJX>VG>M2)_wc%KdCS5M&);<`hH-b;hwfI9Bplutra+qF;YkhM!Lvo+Cnc5qs$;)1n zKzowZ#X`OgP7hqLh%N?+=YYF3kux^LIqn>>G2QI545WMAd$WA(jikqbgLU0fgd@rU z2IwIjNuS7+AoIDtc_%dkULQ&M(yR2JoLmp5k#9{joKc-EO6^)zJY&0L!FrWQ{7AM< z#uPCZ@ftpI$GP^McH=Xa)HZ%K15<6#LQrcf{+%_7{#3QuM1j~n*~S#WQ78`2hMiuz zhqTWS&T7JE$~ z|A{4Fhc#bQZ`8XpP)lH0ihO#}h!TQM+{UfrtwbCBEx}tUcSGIg01Hh0Vo#fGaV<|6hmMOVeMF z0pBC3sx*COaKn8+c^M?hDHZ)xtE8B96Z%SFvo`C~_f>wQE=FjQBTCLjSj+{aYX1LI zV)AhfGM9&$Y8v!l)Icu9b6}=LB?AZP97MykSozDu?9wi@wrjO54GpjgaH1cf!gZF*qQAU!^mf%$YjmIcgS*`pKp zve_%p72|Z=IO~(c8c`n>1a*!6cksz(24^}{nY$JRsHu-; zAaCoZwklO*>Fho2t1cE*}FGWxS5p2Fn z9GTPfFJ*xD?@ycP%1+!^zv}Q{`ZkgoJP~TmCfi)H#em(kS!w+Xvd|xGCWp!Vnlyj^ zql23vnDCi#>v$rr;GR7{1F``p+feR;PmBR>N7v(gc1P9kAanxFd_JHPUr4HPCb)YB zX%ZyKU{R(;&e5=0?MX-j!82U;opUI@F<=-i zv8u`sf^~$uA_7bokQb*FG|m%R5aV0A>hY6F8E9Oaovm8POa2JRazg~eU{-7CBWDK=1eBy*j(wJn4BmjMRnqV9_ z0LtXK%305qvz9aGmPxjo{M`)14&PaWXeC`^C@V?WEZ$S&3rnvs)V}hg2f_qDEu>9c zrA8TAl`k`xhi7WAqz{^Wuz+wOp#Z}kvC%-|LaG}wT^T)laR#C|UhHUZ-nZO|R6|e` zKI#foloVw65SNQia8jE@1897^I}4;35r_%I;Omi4VByCWFt=jsIri1&&^1u;Ao_`>nt%^MlDlVN%}S0rU7oIRY|%5hT_2l zQlJuxYPvReN%$u5^dYRlZJaR8csv8&pH!}oRF3t9`?r;!Hyhu?*bt7!Y7MOfpL z#H;|161KWb{IEYz`i88Lu?DG+nwSSqo(MV*_F%cE@6EtVH)3ggVhUqo=h`VbaU{yz zCj^o6J~sl@zST3IMgU_dX{d@*#9$T_xn~q$D#dx7Slg*SpL#?)Y*Y=*yfIxAv2IuG z*e3ps3>;O+euaY;Plbo8*AY!HkNdPf9`r7)wU5f@R|RW2^;oPI zb|FaBN^urDZcfab*#YC)M-jH{@eJ54@^K5Nn!bm(d_^GyC9cM~aJk&^PiDj&Pa)xV z2l1TzlMSjkmRWhrgjI|smbX{e&=qP8OML4n4{V8Ny1F%@Z(7-!Z^3B>-jnK^UQ)@< z3s74ajHFv~Iu2>(a8?i;!1VhUAy>*w zCqZH*G69-=UIvh1`};hr&$Xn79+WnPN`<8X^ME?DLIgvmUa|rfKI;i-ZK4{?%(N<` zfu-zB&Ot*j5Rwu2+RFXleSH2vf8lFlO^n%8bDo!hZ%lWMiDt&&Q76X3OMuXc!}$@m z)+G&)Ge%w;rW`KIRxFWTIvNoPI|)C`Q>|N25)70C-Gx~F~##O*+)Bob>2G=b5GRRZR4 zEGCYmIz`Ckv%nO2PyXW>h#F%GRVTGr&Lq?a{3{KpLMzKsS}eGeO4W{977KJm02uJ0 zPO-486-7h25xr?ZiB*Te)Q#;h0r&N_8A#jnT1T_>mFtN|45}rt^%hH|7sZ8%sL;3m z%%>8SU&Irgu;+C>2O^w+DgqjCg~np-3Bk^Zn*>aHo4^t{0-&o9VDqCRO_|r1|4Ih@ zs4U8FUtDEW>qotgVHph8Y?sD$df?MOJ4sKeV)w~o_v}c8kxjJneNK_ zL0boN#+y?14iP_TGPCKAW`KA0tr1vVU2__Bgh=_DRX0F{Rw-Ec?jC;yC)T+t!B~6* zxe^-BMk#VKtcN5rMxxXM*}6aVBw7M;L}AOo!2!m&J$@hqCHrOafnX+=d`F$Qg+VdX zD*?XV;pnNf(F^HB3$1JDSbOzDa)km=Ml>^aI$ESbGkSXstT z9q|IFMFJ<{k!1l`*=yd{XW%UxcNOk7mFs=F;*V9F;ccIf5${ zLaJEYpCO>Do?oaojYWJyq()LObo~=%S@5>iE2i?SATqQqUvhQuRCS0b(c{6YrG zKBZ#$&e)?8IX(#<1&R;dnagSb2?kR-l_AP@UMx9+Iq_%kWRQzr4X9Ce=Jaf`Ip3WD z-@32w0x4^bSA(YFN2y|_9LqJh%mQJ|G-`<_4oF2}!1qupMUyEhES$dJI$oDe&ZA9K zr-#nKVFqfwloFzemE2&0N)i&)#h_EX7}zI~QL82t!(<_&Oc{>x%6LGgqxTFD$L^Yr(Oj$<#TOm}9Q5t7Isa;JoIfT8Y?hXOwtOsmHq z&%nMQ-n_80KupdjqKJ)p100izAm5!$Vk03E8HOfGPC+Pp=>r~Am1+JRyh*JYoq<1Q zAeZlKilthI`3_B;;M}a0pB|YHT5Frxcp)$YHY!9d3tH_-7<4_d&xZTb8F)noxJcUP1J|nI8tTvK2{-6gdUi2`F2iH&%u(1W$PK$- zdh#F3z^9%I$7%LbmGI=Ok9LSJevS>Sq&K^0Wukfieg>@eZ2X25b1d8oXaM@(U;5x1 z7Su8PP6kpVB=R%ydjQl7Ow)y)3;Gn(;WMGN*MBMlUTNw0op`p^1#H&-vqcv zimH3Tv(dzBGO(%r4f_c(>@4F)T_Te=h3x9ej}CV8-4D8C1%{=iwd}4~jpmzfr z<)6viTt>q?1AoW>%I**9&r`K1PTJ*p4%<>^fHTk)*M-m-xRL?8M_<721unt`tWAKkbP`zaaN)PMhS|5LKtMt&*-@}q>V{h!*j4)#SE*gP6u83qbOt&Doq^84n={ZqUwd=H+uhe@pnKtKx26N@40HxM1D%1+Kxd#c&>83qbOt&D zoq^6kXP`6C8R!gj20883qbOt&Doq^6kXP`6C8R!gj208wrT_o{ literal 0 HcmV?d00001 diff --git a/api/app/tests/utils/test_geospatial.py b/api/app/tests/utils/test_geospatial.py new file mode 100644 index 000000000..dc82ec3a4 --- /dev/null +++ b/api/app/tests/utils/test_geospatial.py @@ -0,0 +1,108 @@ +import os +import pytest +from osgeo import gdal +import numpy as np + +from app.utils.geospatial import raster_mul, warp_to_match_extent + +fixture_path = os.path.join(os.path.dirname(__file__), "snow_masked_hfi20240810.tif") + + +def get_test_tpi_raster(hfi_ds: gdal.Dataset, fill_value: int): + # Get raster dimensions + x_size = hfi_ds.RasterXSize + y_size = hfi_ds.RasterYSize + + # Get the geotransform and projection from the first raster + geotransform = hfi_ds.GetGeoTransform() + projection = hfi_ds.GetProjection() + + # Create the output raster + driver = gdal.GetDriverByName("MEM") + out_ds: gdal.Dataset = driver.Create("memory", x_size, y_size, 1, gdal.GDT_Byte) + + # Set the geotransform and projection + out_ds.SetGeoTransform(geotransform) + out_ds.SetProjection(projection) + + filler_data = hfi_ds.GetRasterBand(1).ReadAsArray() + tpi_data = np.full_like(filler_data, fill_value) + + # Write the modified data to the new raster + out_band = out_ds.GetRasterBand(1) + out_band.SetNoDataValue(0) + out_band.WriteArray(tpi_data) + return out_ds + + +def get_tpi_raster_wrong_shape(): + driver = gdal.GetDriverByName("MEM") + out_ds: gdal.Dataset = driver.Create("memory", 1, 1, 1, gdal.GDT_Byte) + out_band = out_ds.GetRasterBand(1) + out_band.SetNoDataValue(0) + out_band.WriteArray(np.array([[1]])) + return out_ds + + +def test_zero_case(): + hfi_ds: gdal.Dataset = gdal.Open(fixture_path, gdal.GA_ReadOnly) + tpi_ds: gdal.Dataset = get_test_tpi_raster(hfi_ds, 0) + + masked_raster = raster_mul(tpi_ds, hfi_ds) + masked_data = masked_raster.GetRasterBand(1).ReadAsArray() + + assert masked_data.shape == hfi_ds.GetRasterBand(1).ReadAsArray().shape + assert np.all(masked_data == 0) == True + + hfi_ds = None + tpi_ds = None + + +def test_identity_case(): + hfi_ds: gdal.Dataset = gdal.Open(fixture_path, gdal.GA_ReadOnly) + tpi_ds: gdal.Dataset = get_test_tpi_raster(hfi_ds, 1) + + masked_raster = raster_mul(tpi_ds, hfi_ds) + masked_data = masked_raster.GetRasterBand(1).ReadAsArray() + hfi_data = hfi_ds.GetRasterBand(1).ReadAsArray() + + # do the simple classification for hfi, pixels >4k are 1 + hfi_data[hfi_data >= 1] = 1 + hfi_data[hfi_data < 1] = 0 + + assert masked_data.shape == hfi_data.shape + assert np.all(masked_data == hfi_data) == True + + hfi_ds = None + tpi_ds = None + + +def test_wrong_dimensions(): + hfi_ds: gdal.Dataset = gdal.Open(fixture_path, gdal.GA_ReadOnly) + tpi_ds: gdal.Dataset = get_tpi_raster_wrong_shape() + + with pytest.raises(ValueError): + raster_mul(tpi_ds, hfi_ds) + + hfi_ds = None + tpi_ds = None + + +@pytest.mark.skip(reason="enable once gdal is updated past version 3.4") +def test_warp_to_match_dimension(): + hfi_ds: gdal.Dataset = gdal.Open(fixture_path, gdal.GA_ReadOnly) + tpi_ds: gdal.Dataset = get_tpi_raster_wrong_shape() + + driver = gdal.GetDriverByName("MEM") + out_dataset: gdal.Dataset = driver.Create("memory", hfi_ds.RasterXSize, hfi_ds.RasterYSize, 1, gdal.GDT_Byte) + + warp_to_match_extent(tpi_ds, hfi_ds, out_dataset) + output_data = out_dataset.GetRasterBand(1).ReadAsArray() + hfi_data = hfi_ds.GetRasterBand(1).ReadAsArray() + + assert hfi_data.shape == output_data.shape + assert hfi_ds.RasterXSize == out_dataset.RasterXSize + assert hfi_ds.RasterYSize == out_dataset.RasterYSize + + hfi_ds = None + tpi_ds = None diff --git a/api/app/utils/geospatial.py b/api/app/utils/geospatial.py index e54288235..ee4124736 100644 --- a/api/app/utils/geospatial.py +++ b/api/app/utils/geospatial.py @@ -1,66 +1,64 @@ -from dataclasses import dataclass import logging -from typing import Any, Optional from osgeo import gdal logger = logging.getLogger(__name__) -def warp_to_match_extent(source_raster: gdal.Dataset, raster_to_match: gdal.Dataset, output_path: str) -> gdal.Dataset: +def warp_to_match_extent(source_ds: gdal.Dataset, ds_to_match: gdal.Dataset, output_path: str) -> gdal.Dataset: """ - Warp the source_raster to match the extent and projection of the other raster. + Warp the source dataset to match the extent and projection of the other dataset. - :param source_raster: the raster to warp - :param raster_to_match: the reference raster to match the source against + :param source_ds: the dataset raster to warp + :param ds_to_match: the reference dataset raster to match the source against :param output_path: output path of the resulting raster :return: warped raster dataset """ - source_geotransform = raster_to_match.GetGeoTransform() + source_geotransform = ds_to_match.GetGeoTransform() x_res = source_geotransform[1] y_res = -source_geotransform[5] minx = source_geotransform[0] maxy = source_geotransform[3] - maxx = minx + source_geotransform[1] * raster_to_match.RasterXSize - miny = maxy + source_geotransform[5] * raster_to_match.RasterYSize + maxx = minx + source_geotransform[1] * ds_to_match.RasterXSize + miny = maxy + source_geotransform[5] * ds_to_match.RasterYSize extent = [minx, miny, maxx, maxy] # Warp to match input option parameters - return gdal.Warp(output_path, source_raster, dstSRS=raster_to_match.GetProjection(), outputBounds=extent, xRes=x_res, yRes=y_res, resampleAlg=gdal.GRA_NearestNeighbour) + return gdal.Warp(output_path, source_ds, dstSRS=ds_to_match.GetProjection(), outputBounds=extent, xRes=x_res, yRes=y_res, resampleAlg=gdal.GRA_NearestNeighbour) -def raster_mul(tpi_raster: gdal.Dataset, hfi_raster: gdal.Dataset, chunk_size=256) -> gdal.Dataset: +def raster_mul(tpi_ds: gdal.Dataset, hfi_ds: gdal.Dataset, chunk_size=256) -> gdal.Dataset: """ Multiply rasters together by reading in chunks of pixels at a time to avoid loading the rasters into memory all at once. - :param tpi_raster: Classified TPI raster to multiply against the classified HFI raster - :param hfi_raster: Classified HFI raster to multiply against the classified TPI raster + :param tpi_ds: Classified TPI dataset raster to multiply against the classified HFI dataset raster + :param hfi_ds: Classified HFI dataset raster to multiply against the classified TPI dataset raster :raises ValueError: Raised if the dimensions of the rasters don't match :return: Multiplied raster result as a raster dataset """ # Get raster dimensions - x_size = tpi_raster.RasterXSize - y_size = tpi_raster.RasterYSize + x_size = tpi_ds.RasterXSize + y_size = tpi_ds.RasterYSize # Check if the dimensions of both rasters match - if x_size != hfi_raster.RasterXSize or y_size != hfi_raster.RasterYSize: + if x_size != hfi_ds.RasterXSize or y_size != hfi_ds.RasterYSize: raise ValueError("The dimensions of the two rasters do not match.") # Get the geotransform and projection from the first raster - geotransform = tpi_raster.GetGeoTransform() - projection = tpi_raster.GetProjection() + geotransform = tpi_ds.GetGeoTransform() + projection = tpi_ds.GetProjection() # Create the output raster driver = gdal.GetDriverByName("MEM") - out_raster: gdal.Dataset = driver.Create("memory", x_size, y_size, 1, gdal.GDT_Byte) + out_ds: gdal.Dataset = driver.Create("memory", x_size, y_size, 1, gdal.GDT_Byte) # Set the geotransform and projection - out_raster.SetGeoTransform(geotransform) - out_raster.SetProjection(projection) + out_ds.SetGeoTransform(geotransform) + out_ds.SetProjection(projection) - tpi_raster_band = tpi_raster.GetRasterBand(1) - hfi_raster_band = hfi_raster.GetRasterBand(1) + tpi_raster_band = tpi_ds.GetRasterBand(1) + hfi_raster_band = hfi_ds.GetRasterBand(1) # Process in chunks for y in range(0, y_size, chunk_size): @@ -80,8 +78,8 @@ def raster_mul(tpi_raster: gdal.Dataset, hfi_raster: gdal.Dataset, chunk_size=25 tpi_chunk *= hfi_chunk # Write the result to the output raster - out_raster.GetRasterBand(1).WriteArray(tpi_chunk, x, y) + out_ds.GetRasterBand(1).WriteArray(tpi_chunk, x, y) tpi_chunk = None hfi_chunk = None - return out_raster + return out_ds