From f8180936415fc4bb0c6f44ed149d16a02f2b21fb Mon Sep 17 00:00:00 2001 From: Kirk Bonney <47759761+kbonney@users.noreply.github.com> Date: Tue, 15 Oct 2024 15:32:40 -0700 Subject: [PATCH 1/5] Update workflows to test documentation (#453) * run doctests on quick_check * update documentation workflow to test documentation building on pull requests. --- .github/workflows/build_deploy_pages.yml | 6 ++++-- .github/workflows/quick_check.yml | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/.github/workflows/build_deploy_pages.yml b/.github/workflows/build_deploy_pages.yml index 47d69eaf7..7c4c627e1 100644 --- a/.github/workflows/build_deploy_pages.yml +++ b/.github/workflows/build_deploy_pages.yml @@ -1,11 +1,12 @@ -# This is a basic workflow to help you get started with Actions - name: docs on: push: branches: [main] + pull_request: + branches: [main] workflow_dispatch: + jobs: build: name: Build the documentation with Sphinx @@ -36,6 +37,7 @@ jobs: deploy: name: Deploy documentation to GitHub Pages needs: build + if: github.event_name == 'push' permissions: contents: read pages: write # to deploy to Pages diff --git a/.github/workflows/quick_check.yml b/.github/workflows/quick_check.yml index ce83b425b..c30c9060c 100644 --- a/.github/workflows/quick_check.yml +++ b/.github/workflows/quick_check.yml @@ -33,7 +33,8 @@ jobs: python -m pip install -e . - name: Run tests and coverage (unittests plus doctests) run: | - coverage run --source=wntr --omit="*/tests/*","*/sim/network_isolation/network_isolation.py","*/sim/aml/evaluator.py" -m pytest -m "not time_consuming" --doctest-modules --doctest-glob="*.rst" wntr + coverage run --source=wntr --omit="*/tests/*","*/sim/network_isolation/network_isolation.py","*/sim/aml/evaluator.py" -m pytest -m "not time_consuming" --doctest-modules --doctest-glob="*.rst" wntr + coverage run --source=wntr --omit="*/tests/*","*/sim/network_isolation/network_isolation.py","*/sim/aml/evaluator.py" --append -m pytest --doctest-glob="*.rst" documentation coverage report --fail-under=70 # coverage run --source=wntr --omit="*/tests/*" --append -m pytest --doctest-glob="*.rst" documentation From 6ca9e8cb97d306fbf369515ab10d5f86c65d237d Mon Sep 17 00:00:00 2001 From: David Hart Date: Mon, 21 Oct 2024 10:35:35 -0600 Subject: [PATCH 2/5] ci: Update quick_check.yml (#454) Update the quick_check to no longer fast-fail. --- .github/workflows/quick_check.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/quick_check.yml b/.github/workflows/quick_check.yml index c30c9060c..e13dc79e8 100644 --- a/.github/workflows/quick_check.yml +++ b/.github/workflows/quick_check.yml @@ -16,6 +16,7 @@ jobs: strategy: matrix: python-version: ['3.9', '3.11'] + fail-fast: false runs-on: ubuntu-latest steps: - uses: actions/checkout@v2 From 9a2a18faccf89ebc7e65645021e962b2b40f3987 Mon Sep 17 00:00:00 2001 From: Kirk Bonney <47759761+kbonney@users.noreply.github.com> Date: Mon, 28 Oct 2024 14:13:34 -0700 Subject: [PATCH 3/5] Addition of raster functionality (#446) * define sample_raster * adding module check in sample_raster * exposing sample_raster to wntr.gis interface * add test for sample_raster * trivial whitespace commit to re-run actions * adding rasterio to optional requirement * added tif to gitignore * updated docs for rasterio * clean up for raster functionality * update test to use net1. add user manual documentation. figure and data added for raster example * adjusting test to produce a smaller tiff file. * small updates to user docs * rename datafile * update test to include tanks * update gis to include tanks, add elevations to model, and mention hazard datasets. Figure filename changed. * cleanup test * adding documentation updates --------- Co-authored-by: kaklise --- .gitignore | 1 + documentation/figures/sample_elevations.png | Bin 0 -> 43458 bytes documentation/gis.rst | 107 +++++++++++++++++++- documentation/installation.rst | 2 + documentation/references.bib | 8 ++ examples/data/Net1_elevation_data.tif | Bin 0 -> 33444 bytes requirements.txt | 1 + wntr/gis/__init__.py | 2 +- wntr/gis/geospatial.py | 68 +++++++++++-- wntr/tests/test_gis.py | 52 +++++++++- 10 files changed, 231 insertions(+), 10 deletions(-) create mode 100644 documentation/figures/sample_elevations.png create mode 100644 examples/data/Net1_elevation_data.tif diff --git a/.gitignore b/.gitignore index 49a569569..810a40e74 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,7 @@ temp* examples/*.inp wntr/tests/*.png +wntr/tests/*.tif documentation/_local documentation/apidoc diff --git a/documentation/figures/sample_elevations.png b/documentation/figures/sample_elevations.png new file mode 100644 index 0000000000000000000000000000000000000000..bc67cd6fd384d08170c7a309c5090e7264062b09 GIT binary patch literal 43458 zcmeFa2UL^k_C6YQ)Uk`GC=qc0K}BQ$0j26l5fu>;kgg!Ti6Qi+j!FwCf(QsGC`j); zfQ2H`L`vu=5rPDWA@q9pOVOG0^PKbh-?i>tcilCtl*k=@dMo!_|I`$*~_*n)<(AWS8o_%6tCJ|x3IRg zFxB5@XLti=YHcNSMD&Q@;e95yw%2iD{QQ>R-*CkGhB5y^ogR0%$hzyNHE%y(Mt79i~V&_R&fqlx_@A!h#1&`u*aSFT{ zq{MbX#*Q4)HfPJU!z%Tex{A08KIUNU3%JtNGAGwuQ?|RtZq!S3-^x?-`zoR9HM<=q zT#~8U>f^Ww%3Me7@3@u4jDkq6uxg`Ka5wPd!KWJ&$Flf4=ET)sJytCKe(?{?#>L<6 zZ(Dg~@i&jat(c{cJ^vH4Yw?f!_W$x~@i)xD-!N#s;OF$Z6~8b3cH+!Gm;7ffG5;75 z=AYAp!TjSkF#nJE%TVFbmEdomHmfCH-tn+29J}KorM&8^PgHuk*!%Z_by-?reQm|K zk*>V6nY0h_alz~0-C^gDi$7sjJ#po#ucL{+?<^(oa+-^3A;nNoVPi+M@t)zyn*t{9 z?DlK_25;Z47Eo=Jp;-7U{C?oMYC=4@Q@pH{bKK8AG_-6s!+!j8VfPnn8XKEfKMI>y zDmhU!*r(cCt-(6)T*@rq+W&?S5v^D-e2hA?N1b^QX^!mAl*RvfMy@v_YG^dOMxME2w zHR6PYta2~htj8bBD8QwbIuF;`N9V`Km@xcdVd@rN+qUqUY(iYJS8VyYLuLefo+`iH z?1vAVZGvkb_+=sFX=R4gvy?{N;q}&~gy%>dazd{)Ut0WL6WEs6_|5Pai}g$I-eT;k z9V#YtR%I>xRpbtesvdl;@$W+e8saonbppiTC&R>-0*RFuD3+A>U-pI3nmU3UA@ zW!tQB6czeyyOMPj`r+@dViw|bsw}mwSjDT~_Q>DgEtzA@KG>Rbz$rX@EuQ!9HyDzB zZovJWM^w|U-&4&;`s>qP3d04TZe6?qy3}k)cvPYWea_3dk79qMwH;o(aC&J`b@HV||wzVe>x3 zSw6NV$tM%7@qxVy8xl8qvhqJI06fd%WYFYR?9EzJBXe_WB-e&#Akkhd-}h z0as34ymIHTVGzAty33Ci>8fCBLgIY?{(W>zjM$^YVa$x3Y9*_-bnU!-Z>YVY@)vkW z#iqo?RI9Wp9`CrZ`Go;-v~v}{s95pd&i%Sfa2MO<|4GI$|D=o<%s(mPKPe-4jDHHQ ze^N%sq5jWH>{O-zu8P9V65n3M=x#R8>v^O#Jm+Es|Bcbz7U=swYWJggMXuAmcxvP5 zAzoe|KQ?i>t|C`*>f6me8~BO8m-(@^HJmjTQc+e;drf;$M%;NsySX9fmk$RP0y_qi z{j0FBaN70HZyt-`%_PQ}4d@ReXSn2icQHjyl<;)77)M(6_XJhlT)X!%%Y{E!C#CKY zcow08D?VdQTGpj2*6!I_JXviW>U0qE;mcAs3t=|fJUNg1?#A5#_KErVVml9C-fCSm zY1HmI-a#KKp!Iep#Kn1K)QIs;jnmj8C##Il6l#m${3XYl<@DOJE!$Qm|9SCvPoPS7 zpvrmEe!~+Fm+*_!mfa#2m-je)x`{8cpP`qR?@-rflAdmOks0|uL&~|#(WsDSbI8_2CN4~G>djm_ zuI}EqcQ^S&`%8({NHD5Q8sk-m3YDYfN?!!)YSO0qaj#aaEm2XD%<9bFv15nLWKVE8 z*dPmoQk%kf(v+lC9$Yl?B6lFYGD@%5jT&xEK*y;l8TQD_+uPgI(_wCsLXJ;J_AZ{E zrKOoQ#uLqx3&&cMxG`?jq@3>b$|FAVdmXtR9oD!=InE|#>&dJnA#kK(5;7CV~v(0yeG)}?=7qmp)xV~u<-lA(~80p6gl+A#p=Jl zwR*IsB}&%kShpu<_nQs;z8&^mu7#r;CSBbs`Ep9)Vk1Ja^Pji4d#ap2Zdf=jo|&Chn#y*L}ZrrP=a!Pt_doX$}7SnIB$G_a%G zk&kD|cd=HU@Om;%$I;R8iCwFyO^ac$aPCemQLPgTuwklfrBQ|>`R*^7-RRCaS8wSslt(G^ z-7Um=$v{QNrgTRdm`HiGXqUG21ePhyq9r|dcC5|d)2o~4vX-`m2`NQl$O&x$Hz@X* z))W;Vf2c-`nz&2w+Ua&YNoVwKF%#1UmLtB@HtyQGMRMdqA%)GRgvn8-e;n-D ziIQ|yB#aPL)>hMDAI5%iAtb&mdm)GFVzR=Ycf&ncDdo zrRb)H?QV07$nsY=SCPlstW!j12$yGlxQnM6Q=;7G=WJ>v$1Q9N>E!bA`@fy)ZZ%Iy z8qatEXDzvSZVH(!uK4r{>OuvL!CdhX0ZSs%v)il+UYn~WUQVk?OV%kUOD-5BM*8YF?YdmYo6S=S^+`JR zgw1SX7S{{!sOjji#Xh-lM%#s6s;r?AGOkOh2{-um={539GPMzsjLw9_MBjQPX-QiGEg)I%~l12{SrB2wMqe8X|TKHpI+ZKpcRh$oxHt9mfDoWsHrf@u&I+}DK!k%X?lqd z`dG@Mx2J1*ngr3g9opgg$`A*Y{rRT~3OwQ4Bc~c-IXJqX z;A$$U2i)MR`gP=8YmXGcjSgo=S#}Z~pE&l0(5tQ8gzR9=i8vmmh&|Er{(e}i+}~I_ zBi!dE(@46`rLX~7)`MBqkqA`fD3^{eC4pPf2R> zOYYnh6-P$r#b&sIN7TX6+nuPE;zd@dc|R`%k}rRXKhHXlx|6@J;YmfdA5o4B2JZ2fh@sx1<~ z&rR04_m&Z6YiA6~{Rm+F{G^+$;9(uUJdoFiEax%}Z{O)??9#)lSD_;c&%8_2O1qFk zniC>;Hp6-7_Tq`~9O@_#HLE5lle(sxSoS8@)`-*gK$LbZq@{FZ5uZ``y5Jccb@FGY z$NkZ+Ln>9=DeSvfQmjULex`i_ELR^z1AK(|^-r%JAG6W7vQ&vZuS1(6T~5=NwWW1R z?&QmgB!gdM33q(myuBDgLB%qUXq1fQf}bWiSI4+TiuUe)g&YK%UGr~J zQ+?6xX{WYG9^DIx7z7d7aMIz+sryi%cVJYD;>z!CW@Ak$9Mdc`kDuN6$S2ma|oXzZGEfyk>x5rwl3*^ z)r0n%U|jW>X?~(yx?gf04rBE!?dGbrBFAQ6(IEMvZgSu783J{xkq@Ut`x~1r1E~+IjSL!ID(BgN zf5L8`lS`c2e-9lqQMXSw*A5pIH)Y$s3#~RHgF83>FolA(_5?Onw!r_k|oKQ=8|US%+U?JN!Dik6B5;j3@{B=ULsix1_lu_o_H6Z6z8 zpymmz_-a$m#!@htt#!Lr=^-+uPOz-@iGhq-s|ynaWvv~Nus=hEsDN{@CX~dCyvOc#%(}ffWg4OqLoAj$oj=~+$SVbo zg;K`MbKaLl&&8$?=RWcIo}<<^S}xg-F!jRYZhRr%c<-&Zm(0{v%Gm7i1m@P1IFYsthFx`36>448(!EQ zYrLNnX-@4Zupe?ShLgIRI#ws)LVC-9Tw&38A#KO=xmLY#5zE#R$Y-#U<*l{D z;lE$~qma_*Rx@}v-vg?jF+k>P2ufmKpRY-Yi;3YzCVY(Y1c%(U7erT?dz;y!5)&=L zBg4V)j?-F4XJ-OGLVfbYejYLxAgKSPDdxZg2x~i&ryrJzs5$wQGm2AalBBd<`}0N zWR$%T^~asZ1}4fxiA>9aI<>6mjNA3cNS%ctDj&cxYW76*4NSBpX+(XZXKK1McD%2mWc z$G~l_HPiH=sIGIE<4eW0i%s;faGpq)8{t+}r9Ld>xi;i(s3hKRIFH<}-b^ZCWfopxp6?E2h@q^dRkA@X=(5Q>iCk_xQVw zL@%60G4M|FwMqy1z>`}O1&Uu@m=0~lHN}Ft5t6gyPKnc6#J(=q{vV)5-zDL2edAt< zV?}eLDdku1vF#+65qvxwJFX*-OPuTtquRc290SK!7F;~j^D;59itTIw7kYx}p#tjx zz$oe<#e`7XBetR%0^Vb4f5&|YI-H$O)TD_DCAV$h6;;A@!p=TA?1c{WUdJyFCSGrG z+bRwJ+a2X|;f+k0F5I}4rVVeO#i+UclM&h9#*`~(gH_DOLA6n6ol zogXOjI%awqG1IJ}-it7d2jBS|tB!5}M?9aMy4QQF@t{gf$xG#EeMs4-niBsP%*wK>daa zGlQzf0uU)@I5pvWUTSWt%g0lZiQ~vXNZTHRC((^HgM@en@{$@tTc#;9c=-{ik)S|85r@_c`P^vkRDVJ=j#5J{ zwCgVMA`}t=Ud|G`;1G#8`L&WV%}G`dn|=v9vkZ}#P7mhzL&zyNjSYV!K31tWCkURI zCIbKzm1tBEu%G_n)RvwyFk;IYI7j{fPKmrUI+bIWltd&-Y=#=)7#IPEKE8Yyu&S+) zfm}!T=NBSm+*{)c6=OG4CxotFWz7CRD6X6wgl=AEPE_+@LJn0+c$GGW6(TS4}LHBu6jl zuAM=aiwx~{Q)Urmr|64T>1k|+;#7McjwH?#7;22#TgZWn0UV_dz_EoIo>ahC_AsFhO}b^sa>bZf}LTs^ zmM|+Fb-TIg1JRL@3*|$%O61g~p=P+%ZYaAF88hd&=cNa&`y5lv#o^wE%S!WkhJ{N> zuJ_pN-c3h10d8XO=FfFZY55Ol!QS5$k_zl;_@>Lj7w2Xy-j6#PKuWw9Fri1-Vz67P z1dD?!w_@mal-^qrJM7wIX`O`|^5skpFR_cXGWp%;KFEo-I{;)um!4-x<0G`hWi$)0 z4zB4f7TWH!quAniYY^VmnLC!9Zk}mc&uvWvtY8gkV!&8BOtk3pU2K!%PJuh>(^NNz zszJ|(4#qwS5Y??8b*6;d_u@0$p*}PJz(^7uqfh`L@`h*?dy=s6!qMuteQ( za8CT2isqp1tHH8gf3;$@7S`g&eR3&umrKvtbf#%pbl)3_Gj54v=k$OfMpg8j#Pmg?p zFDXN-M>Obyu(&>7q4}rKNyCs_m8%1|;V_s*Vs%lJm0c(j^K7^udVy1GavW+M-_7Jm5-tJn z6p>Rn)Tpk<25F~8bar#Nl&jMdsAK3K&+-e$k2S@kB*9n(tqrQlij)~d@bnAytx?tE zx~Gh`EOQnF63i;wbJu0L(Njt1I*%snJTw?+W_TaltBfi)IWGhkBn9x;XRu| z-7UkWps3WF?N>K|a zcLGG&`AO5v%&fcCeGZ4fkTgn!bD0us8p;8w8y`vmOpyCECQKit;1E4}{MZ49qS|UD z+w?uV*Gu3po*6zi+Mt5U2A^Vp5(Hp1+NnY6x*C)zCQi8$_$|bKAmHo~FGzZYAcOR~ zoO(VQP@LyC%>jUP3opG|*+%y+@bWBy9XR_=V+yhywYf8SUXz^}{_tzyRKHbK<@u??i?m)pg62eScm)mqH!nVeX=Ut<-2O(5GWRzvOhddjb zrXsXd2J)F3J z1Z=uXJUuVx-7uUPiIvtd_Gy9)Qv#*G3t+iE{@3XwYggVHu$5*AI4DdDm^URp1{)j- zO;1nHg-%I#Ti2XjRhVeXYjfnW76K~!+_!tmfD!}X?2XzZfCa)k_yItQg2n_b6v1&R zZleiK6zzJ$32>7#P;E{=(z1#ZauNIFv8o%OR!z#p7Wav_k0qV$v>|`WGeG4Sf>s|? zT#K@&b>^iye|~cZ-Abv1`6smu67j3^5En4Zu=6A=Rr5<>ZE2FuM7c@f>i`7op3Jwo!QoaO|2{Q8k zH_FQh*x-U3l-}d<_)C{Al|U}cNUpvnNQXrfS+iv?QNi5fkGd5LT+$PVt_g ze<HIDXxR&_#)98U*5^h@7{~660qqBQimpkV>tqt(Do9+w>y5OL7B!5 z8`JV>+*W5j9XwF5o|_{w0>`!z`NX0elAgO{rwgMR|E9|_yM=owpPG!$`q}LB&0ClJ zPaHic#brOB)c3LcdPNd)2=sS^xO$lgp%U#%gSa{>gCT*s;fobbR@MnXjRZjc_~Zzz zJ2>ThJ)vN*_!j#5_S$lO_sL4Au0>f#B%1b0I+dWz2|=<5O7cR`Dk|1AGRiDx96u`X zz?#XTRAp0q8<@$&aB%S+FaUNkU~7rv$4e8tCMl8Zus91Kg46g=;cEuO^}P$d3h0mt z;CicwltvL5*EWer+m8x_P-6o7Y@a+z>`DRs?w2$Uv&{BF{Wi1)&6YVGBJt;+71d(N z&Rx^G=N|w5+j+2dyqm%`OYoz9sNsjw0==n^!QR<>n|CQ>3n<}%Kjp!93GQ@+N=m62 zhjMEhlyH25w=f>nWpGZqAp!FGS|3}UZdj?0z7Msf5GHzmGYP@>6Lhk{$L9*@_+3uF zrog4Jt<54s5p|vr3C4fk5&oHj_b7b_&f)St%a~c~!$nn_`T9>ro z?n1*mBBmdoJIs!f3~HVnPrLH^mPu2hMy7R_6s!&#f|!kF?HkuIbC>w=7wtmV0#oaQ zToGP5k_Pc3rCy=nMH~cQSno2_D}&1*FQJuVg&K3E;3H=br?xCI(o`4OUc%oJ5UUUU z{Qr)KT^;JQV)4=LWa^iCKbCO3p~BCp1g;|*3QNVpjaPX0Ef$rSi-vGy)1X`s9u*uW zjy7otn0$B|qq`<`N9q!W|2Klso1GsPE_@>1RaEF$HbHHF(xmcXf4wH6^2bmN(^5*b zuLoYK6Ef9mjcWP$Xw~JnOP_+Jezt0n2!A?@i}}Gpo3q*%nVE!NB(2-J`~qg5uqXM;7fEYV?{2=b_)LGJdGI}5cv6o`m4UK>LAm>rOT^Ew`qH=bbbn*k zH?L$}xWMVGa{)3jCC-Az>0=&iLjQ8z4~vc$^IVfVV@ma7Xv@c@RhLgK!3zIo9rlNC z2^uPxnj4+f+Hz$R+pEQU{F_ZJ_4Cm4B)SFVbDvQ8Mqj9J05IPswkRP$53v8i!Lw(y zgtS_nFy{a53AQP#bwA4WWqW#)J8Ef3|3-Ycx`guXYZMIeJm)ojt;Xc_s!&gV_T@GF z@vH$}?_lMO@TSxii9dGXFC>b;UW-dG&`1@R(fF7i7rfr$ zSvrnKZd${D{LE*I&&|SC`GjGkwD@@LL%oN+)Ol(NY8%Ssvt!h`cF&(%!Q`xwcHG|j8~as@P}xz_HgbGYiY z&P7tKsmkP4i(y@S%jPsi#?`vR^+kQ=tFm2MEEHu*>mRtfpDXFsz!p6v>d17>a<&{n&bG76rLOK<7oQ3sk7) z(fM4*ecbq3gKRgcs2%ZBITTSr=mG#}*HM!T3M7n0_>7qlCbc5oES?;)rOTHh9h&<3 zlRd%W8o5KZo?Iv^q))RdD4A<>44`Z=D%cHwCrb>8*7x7sfuj)=5vTeo;ZEsnzVzIeqY4{i-NYhDJ!xy z>)9OzBBp!i%hGq_urAea=qd&SyGdnV+Sq-VM8jDI?FUp#cC1XzKF95KPQ!G5Q%J>x zDDQ`(Gx5J+QWYMltCkHu-0PplPO!c{QIHKcJlR4fXOEYa49CCcNW5RIs9N3qIAf+$ z?r!n?148^n?(pQ};^De?MhWYQV*56wXLR`uJPZ~-GNaQrJo&(NuJv8f^!2MOYo3?o z^3HGfcb$_btam4k)n!e?tih1?+UP?CpZT_lUXvYdfbULg)FQ5&(Xy)%f8hFAovD#m z)HRvzP8e|Z-+o2)M5ym{q;zQlOBb<%afwf3uR+lSdE(nw(e^VxE^peypcUwqg5F zuvfbEyl^^B)T(^tK35~jNll^?oqfHC-so@;y}>59)Y~y^uF?jRsE547M9}Fg9^BpJ zwV37u52Yh^lY149Oc1TB+@-nsLN&fQnmSd+bCwW&9B9T$B)<=^Z zPMue4m|T%<7jo_+7r#f$2gFo<^(ZbXDk?7c1hbRd*XL0r2dh@>1IrA9I_d>h!@+Xd z7~oW+;wYmQ5B;4hL{t>2XpuJqP)l;zs-sC3$mQvp0PePMDyE>xGZ%} z$KXAiQkhiooK;}ffquPEAhE9pfmky_w9NbWfB3U#nv@W&RU9FkL{IS|cPGfRt2M9& zLQYMuNi>gEbgYbOax@y~>Mn_PQ*B^<)oV%%FDS!FACKz8&y9CybcarN9(H28!eC7H zG0XGJ`e;U*wRk_CFPpBz*WQ#(T9umh{u+m%;fs`YAGW*9U@(M>o$<1^Uki2WsJ27( zlJB)xM1OyNbj?7!WNRkhTx1WAR7*jhut#3o?C%Z(!SBHN8|17TF(H zpMf*&$Am?muZ!ZsKRO?&f_VPV!{+~3gE@NU;G3?Pm$LN0pl`Kli5lEd>DPKZs%12| z;p>g&+}CAmh8(D-oXd{S=fY#Mj)Y%YA%;XZIsPNzBYG^Yj|Hvi1d4I;izmrAaNy^w!8UE zP?Gj8fzRr)DU3$i%mKfAhwF7zSRD;}Ni0qK+r22#S3!%|C=#Yw(@>BpB>f}l7ZarP z5IF-NahTJ(qt7~|H21PsyCl67UW0G@ZOW=G=cVw0@H@P*qM(kie!^k zSQMGXjC)?zE@4|6klvS1oqLcuaDmAB$|#{ty{)XoXYO)?@rd1wOtDFQf_YuFPt&vEnUlx3L~nJZIjxssv{L-?4lnm-y5Pt9GOJ^lO} zOSXQJv#nRBUOU{%nh1JnRLz#utjcQ~?*$bU^u#=F#vrd_lBrg$S@Nkye7ofwRPiiH zVOa-SX99H@=VDEivt8?Yt=lYePH%J1x|7s8HZ2#J!~5}z0@td|ydT?NJ<|WM#Sa2e zR8>d1*1f*@tW!RMqQFpLpAXsy?riPou5E*%IpEGL!*aZp#&>BVcgFU9+=nZa#ig@1 zq^2%A&-lihrX!}EhXQ`nL&6ps4%@xcW>wE|XY z`8}sZFcpp+8hCD5ATQWgc_*B?S%2{OlvY}`pX3HH0MhmO$UR|dAJNt+4;YXoO878Fb`G?HIA*9V!79D z?-&@wGRL$VlgRr}fkQ>DTRrAjt>4<_`$md|K=Y~;O+7{HYl%54AvpNLfkxg<$hUhX zBH8_I@E+uaciba>4@=3ew)fd-Sg8N`eQZ>EZ?v3FTYk&=x+0gZcU)q_9{Kx2j2;ZQ znO$=p*xO>0BYRly!#AJG_nf;Ou9bM4;#8QaIeH-{yUpIGR&v_MVcOWf)I4}jAU-(7 z;#-zwHf`th+?aD6_1j8o4aj%rkyl$ErNd~lv6x@g{nn5{4$7mI{&}Zf-*7Zutk&{@ z>E2c1V}4h4?0g3`?GMquH&slT1Q-UIT$VoZ{9OxgG$wPL)qge$wMM=ajGKbE2h}`g#j5 zK|@2j7sKD5O-Wp$vtW`YFlX=XV@PwqOvzKw;LayB^R0=nkr|P6ntyk6Jn{|H$Og&GZ>~t=rWV7$x0BPZtV1CaY)gUIOx_AvoTkA#vv1}6-1$``+nz?clr)J=& zr$28@iLvy2Tr)mLw=UT=;Uy4^^DmYn2fmZ9I-s!+IHT7hL$RCw8gJfn38_S?y3$GnD$6$nD z8Fe<(H~n;!pqakWVS+nx^5O`4u8Io(jw=a|AVlr&{8V~Z@PrC#{mX(*%GNAT5&cu) zzr!cDjk&zm*Nbc{GWv%{;CrLc@(s+C68Lv?^rAT|u*~wCEc;8H6+k?Q|{YMne{4c6NE+&NlK1T>ZPych} z7Esl%9*3$VcIoA+a=AW(^0`8e`{Do1tNcZH{-ZX%mj2*c`hFDXm+Nm36D(DE`x*nd zI8e175SYF973*IXV6hQ~3fXLrMK$|HzK0sKhwv!Ag3Dg`*|&fBl3yQP<5c(`-|~ZC zzr^;23OwL{USx@<0w}_diZp-=tH7-_3`OV3A9eb2@qCf*Ck6KsMgeM(WrX9tS;#pN z!4Ql8T#=$D`T9tP>xGrq*)!T(K6NcW)4wd!|Jma$<1kMT08y0*5!ax%!-J^%z)Pu1 zP)}K4paO7tz$B^SAx}AUj8hQ&eo0MD4#-A`zMBT9v>Xr?HmJo0j{kx^wX`KC5LONp z5&JJ}Q~$Eqe zXv!>{mSfo_avn4RxnBa#MghwKtAoBH?9Cbw=)up4`27kZ0q4B*6td&23zQKv4KP^( zV6jT2h;EDSmb&Q-Cx@NClp zlmmsD6(PI8VxMTrDGIV5gys@6fXbG)%?ZE}MEp3S0Wj_4&V`BabwZy2bpZ9}W+oBC zSqA#GILcwI?8AP*$hZtm2qSx?3kS45EW;^xM$p0%L_2peF#+4|NZywKLWDWAyb)&H zo(p7p0lR@_)JWq}2oUcC>zNp|cE=(bSLZ0|AHtnqxs3;68bJ7njDZ@yKql~SPdmrW zjGBa?$q<vH+7j!?X`y-y-Up{$S&w46@%VBdO+F~XWGitoy@o>?-p+d&z zXZqEP-;0kSS}6bhq(Q`Gmv1IB#jdwGyl4>ybS94C>~E!-?=dv|0!)<|rs+?@#%$P$Qo1T6Xi>TY#v} z>Z(Xe*1@t}Z$bocI+=a`m}t?rHye;2ULE2L8|O%dUiZ=dDK3R5*k2?HaiSzrFC%g* z+ZB?omJYsOOgGW_mj3c^m?c!Wx7AN#_?=v>O~w9T#AjNc9%1dq2t)ho>ACCJ^(!-d z^3sRi5t*oK!mt_D^s--GVYQ51aXa%LVIIdUlX^Y27oes%cq?Fy2#um9_VDoVt?FWE z13;s&ssmcMMEY<~Wmw|*IP;Uwmbhh@i+lF^WYY;Bl%gmC+X77M%Rm7k!KAqx-}mjbJH!dNRSD`sK0mORyod_#_ab6?Fx zvVfb#KVA%KgW7UvlWPTvfq)9gK#@{~$|U!ZI(cPP)gY{~U`QY8Q&R9>Kkf&inuXoy zSZjh>(t;2S0MoWlnnFfZ@37#>UCVv&N*6YUeaE=(WE)W_o}TyiIvoh$YWYNcdVBNx zFoXVvIO-DyaE&S98E^{}LabsK2rinWis7sAQv5k!9Cgt4yVb-BdmBy(0yn#B+09yH zw*J-+9X2d}7N`hK;(96#hyHT1J#2Z;hW2I!A++|Wb*F@AGN{dfzA8L11uTgp))KH> zPPIw1ouCE7Y8yi@4Jc6$u%hr*BzahJ6OL;Z>5?bbcYyOm9zjk=$f`Z}7`SvC8Fdw{ zVCRWA)I$}u0;?Z@;6(=_90f%u0~_dcTBhp44nPn63!6Adrhp*6L2u#@L?x&ul>6>x z1`Isl+S{Ach&D~6pIpnf9SEXOv|Q}`;cp?6Rp6)y%;wDb?dG*KHlP9#iefY}X8K;j z1|zp;@hu+_Zyc|FRTS-Lc$zwbL?1x<`?5fC8vsu*8c~c=Am9<|!)az>qR&m!%8~r6 zCc}q%1H?@|1bM$bd7X^2uyZEMdFqvi(3Yjwfj?=VfQaUY`mco~;9MXQAj2=*sXy*P zMFYE@pWs+LA4P{iwg7|f8z*8?uBJR%m8X(%q4Jiy#-J{z>kOGN1#Vf zq|8sH5Dj7@-JC!!n9*sToZlbsF!ePSM;;q&D+jiB&TCC14eD&FMpQ&b_MK!)5Nz4N zTfp@LZ8;h;E^3m>dk7#4y-T?CSG7l`p(M3lB+)~ps4|ppfWSl)Nnm>Ledj?GS#yl+ zyTBVrE(Z2M6u3}qVm&YlC+hs&0__J-XcnrK5aNH51mzl#h4{yy{|}1AcWg9B1`usO ztWomhZ_5UP2}^`a*gzIkW(+i9Bx8_~ITtP_jFd^JX^!+=Z|-a`fb2+FMcaR{S=&T=T%ty7h+C z=~3o|wZwQY=-S|TnS>{k;N4|Q+yV`r1GDedi?$(62baRfPhv=afS3Tl%3=5+o)IX6 z*g+7{D)!kLR9a-99r*$C$Kny}6bIDm=;tG!m*=^Qtt7Q4)2=`W^Fs!&?M%4Ah?S|uUnJmXxWURE4p*Y7GJ!)8 zPq!`?GbdzZka*UyHoN2U0x=VqLSBejjuo%=iRprb+WvY+f&V2bM1kH^S}pyhKYdIRC%Zj34nDYy7Pr>{UizadIOd_AOUc3CCj%_EjBR<%Z}B}ntMl(lxpUhTVlvtC zUL^mw-m`5b7n9-s&VZJGx%;0~ny_*lD%}yspzKzZ`JV+VT8E2j-G16^ul}|UKfTiP z5iP;Cm8>Ql|F7&&(qULv8zt=!zNa^1mSG5JIWMP2Umq!)>3R?bn|JiTb-}L%oIXR` zkBj@S{pJoq1sM&moNL*ukA9F@{-VE>+Z}?+Q-0cf4kwuDKN{*kS*ry@Ow80e!SkZ? z2J7D{_n%tOQXBhK?DsbIGNJ$*wAjLD{?Ddw1$fXwSS|exCQnV5KMZVyrG7W0CxS0% zm?L1$?u@VBxd1=8b*qqTqX`YraQFDnwtNBv*HPYr~don zADsPP9iRU%ex>QH&7UrkR!xXTN-ai@e?!DT5L_S;m>7TA#-(HAVf$|$Y&kpo3;fg} zsF=Fi77TeNOzAtbVmYw~L;db*D!CH$UL>#k@XODv7_c?)3Ax}8J{pwNkMO2LY zKO=`7f+0)JTgr858Gb;ORtguV{Mxw=T7vD>FA-jT?!SL#d&}#Ux&!g85LN}GPJPH) zP^EJ@%S_cGqgtdc4@W^l6~xa1e!ALyCy+tVyM6e}fMve}n1eom!62>jg8CE)N4FX- zLu}iAgx3>(c>6#=PKAEs6KpfrS#%&eU6d`ofv_w1!&CU!&arIc$7|0|l=7P}FjdjJ^c0zf|p zrfaPjTLP+LkSOrH^z6u*T5pHz5Kw|xnUM52Lq)ouRtWiI8B#C;{mweNKX5qWii+jj@~N{(M>I)M zsOWRi{R2@D)v-toI|KA&Z=iD?M`i;hx6Rv4qB1GEg~Q2(0923><9CfY8crclzK0os zdu&uNh<4zn2)NJ9AjR$iy&ZD7GFcX#?IBs&@7>U}2*`)6zs-WEz}*XtC!jS0+nf!F z3PCJZ3IgGm`T8A>9JPeGtR`fp|6?Ct>#)QjF#^=%gN+|myoXAJs5oS^`olE@N^i$V>Yr;r@)S2)ser$?4%<=VO=E%2N;cb$fe){ zvUNg2a{eJ0c9h74!& z>2)OS8x{MJr@xi*;+fy^-p3O|_QeISzlrL1T)_<{`D3AQQUjo?$bgX4cwY?9#@_^{ zxm*+9w_0coTnA#f-h!+Pfg{*~E4ok|cjgj_#JW)^7axe~7P(BWp<0|vy4w&3XdqO! zmX-6%u|pY+JOSIOE`%y>Cz5R;>4;6MX>1QZS&C65v>)< zBb1zy1i>y57-dRIU56k%q+daz9XiZ$kdY|?4L#E3A+R_Z=Qel+2T~eBvc`01{+o{= z29*VM2smVlj=K7R!h=ij=pPx&^kbk9yf};0g^uNDB9~b0An3O`9i~RU5|Ty1tB=;? z$qm-@l{JM57|7j5k?2RiqTwqN9djw%AcMZK3<3%1xHaGyHCUxEktAh2f`)TSZ@~N= zIQpCM5e0EsW|7h`w;`&*wCcwJhm0^_-eQ$bx_bI&yyXfp+rF(J2Zm+`M4@Dy?%-BD zaHwG#Ov@)R@PzNM<3MW<6;7B|h{-Z1uuB)c#`njo>S471isO8^QXM9$M`@nVv8O28 z1sN83bm0!EhVcCHE7zYJ|2%m0(#=P|tyy8J{A-e6wus3d&B3D0Ntvw#kJf~fmm17{ z_#S@bSm~|3x*&woT(ytuYL(2J`yo}TL9NqsOc$LFsi@wlk*Wx$*n~@uC0o1G=Cn08 z3|CG^Ah&@az%goXZ!ePbix*PB_g%7ah0cPttz~(nVuVCqJsvz0v%kdf=~%NJj@&P~ zm-}JaZwlPIrns(v$3SC@yq!y=q2KcNk#kl6lZcqr18d0^czh1rLSwBeiUPJ7k z=Z-`g)6!9tZa_z^^NSET8Xstu6wT&qA3+0EGAfL{C73VyEG6uxwc_?DADPJ5MDUOT zER_tI$lymIpPhW(AUZydrh|<+Z`!npeuv-XYIH4{2r}2BH}|5u6HUc}k#PNTk&}?e z>_(Gb4w-4g`@>v|m0Y@+(Tt%?y<%uLEPl2l5AZ=k@pM74&!b0=#+;FQ5gK@Q5Ao=B zDzh1Fi*=w!d;#=V7Ut7ISdv)CSOFg0k(87~O{!!_poC#Ty#|n@Qaj+1JEZ5|uoP`w zyXX2VZga)ZqX$8^L~<*Hgya}ZAt8sFCwmE-H)VF|L2qwCi}=9b<$3p{l54+8*v$D( zJoHO+Iv2ed<57Mfh-oRj9 zrlGS6iG#>Q2GWlJnCRpB+rHy&wV13g_8VVi8$#zLZGLV#M^75pH7Q`#^{lk$EiR>O z1#JpPqGZ%cX)$r^{0f1{51`llZe0}>m5o@D!^kngEHgCt$^kr)J{U2!oebkP{2-8+ zMUx}vKZ4q|8#xS+xv93(APGddG(yHGQO!cHn*^5>5c{seNGlyQ!^#F? zvMhY9NI}u1L!}@x-W9BULjh@|=wFrG^Xos;#MZLeUI2Nf5WO+QO&A6f<<{ZR6fLB1 zvH>xAIe?YX!D9E9SdlmftIp*Ud6d zn`RXQ(sDs9tro$h~E43!@6*fXut`Yeuqj3 z6H7}pEKUxjECyfRJ)mdj6eg*hv@y=g&MrfVJNib1rn0iKE>ul9g@cx1v*$l2qZai! zm{J7`!v?3=4=z%s>pZvvn%)g|&fI;6I2VN^7`v@KH;%(lt;lAm9M>nbG&k=;jqyGB z24&%tEbu_qHutR6w>$pjBbaT=co-ag;+Tr9qKDYGrL%|G$Fz(qWCoqV_hF9g00@Tf zJ7l-*2N&m0rX_7n>zN+U8+DvdM?b9?{hYImL2-P?AH05FRb@kIM9;-&(v3VElEB zFWm!_Cv6nd`9u1$EIq<_$n0927>RkFbjj}ki{whpEqgDpyG<`1JPg-&&Q|utxCy`Z zi*e0P80#OH^RG+64s(LTcZtyi&qliXOXgmQr&2Rxi*-Nd^bhQr32t=Z;-*3?S(|T8 z3@+n&zb>PFt!G|{*BPC1v<>jjX_%2h&y$>jk|7Oz_Dmi^5 z(%-=HF1N;ENq(~7nKMcA+d-+aciCECWxkWSeqA~Z9?yFmcREJw(Ea=Smiqvi5*)+t zHl(b(u;oA7^YGze30^?kB$urPmT$=#_nA&=%^w6!)L&owj(mMzS_Lg}jM%?k;wNvb zWE>tQeq#EtnzZV&{C9l)`-xDf?3z-l`}=A6*|CZ+en8Rqmj(XG`*im@7V52dy7kA? ze|H3$do{j1_}6v&+4??^kD%OvEnoXzw*1-yc{K+`52s}S4s-vzd3_xAe8*PCL3|~8vVEa> zLVpcV)fUDBZ4zzY%E!k?oh3qSFX$Xi&P7ZWpzUaB^VBlX;2APtvVK5ex_aAu%_i3E zNXt#;)tfvsGSuPhW4;5=-ORe|+_?HPB4;fI07Ias`MVxQ^cqDp=?#)d8z5UP%uGPZ z6ySnZ${h$i8-rm^kixp9y?%OB-|{-7{^Mn6vI6>oFODP?#{RR9N)SW~#eS;^nn?<^ zD5A~nKuy$hA%ZU|@f_ab%ocphe)MDd0a-2GL&VZ$1cmf>N>mYbY89^8`jeJ}}fMokf`%sDcduffRsIWnsbqkLO7mfMq{HgIr=P$SGrC~5A2}w>s^2yV0V6v; zTjy!}A-an$u=*SjD2l*4ghnSr4TUH~)A@*`h$v^4mX>YpAx~MCPr^DZkOmiqBLP#V z%zdWZy#$6s?E*e!(&c(sWT@|(5Un?Z(nXVMK=8|j#yfy(ou+*&T1IgEB>}cQmS&Bn9B9u* zuakifX;K+183MA8Gc;6|8rU*W(&G+bYV7B3A&@G7;@}cY`L__EKOxVcyaRl!Pld=V zNU3xJA+i_4_1g^@T)=@Vo0Y+ARXVTm)~)qqD6TIjwln>`*dMZ2=YJ?q|B)X{C^m*<2_AnvOpGcch9y{2n=COS8IbGHY$Ki_EEd&A7yehOW~ zh#iyar3Gl{4J%7pJLZAJzRbCM^y@JQb$3FVll$t>{{8!J!qD)sfX|>-$~)C*9@q3e zAdlXfpr2`ZE+Oao$}0t+S3Z{1YAH}1W*aBk#c6nhG{PvI!2(UjZt(o5Cd)x7JvvEA zPaayBwjQV$IjS$;lJ4KJBp>z%cy;t_4xa=Ym*ci_fd8r2p-)XDNFs+C0_0ztBm60# zPXR7M%U2@|v;|;I>6-5FsK?G7@`R*pC)letBR6hWeOC9~EFczf-;f+IECW{Iy+1a< zP>iC3L#81tb)r9jB2da132ukZ2$N z1ynfnHi6ydPvCi=P2_7p;{_Dg)r_IYO>Hu+tMw>XLBJktLhI}6 zy^~v$>j%J74D+_oc+J`y5p(8-1dZ&0_=N4r;002mmJA$ApGfjW9^v-6-v+CLD7Tf|-iCwNHCO3%nFD6<8 zYX`FX%7z;T4jr}%71qG+HDF{V_-MpSbSIy|=uU@qx~*#gxc5gtmxDI7dl4AG%$bjc zLThyAlA3G(kIdi6zHKr9=bs_9Uwppfy!#F&DmM`|?Cjx(! z`{pqkOCYgO*I{aMSJVf|6uhp+-O#NRKSGEH)#OAo+*CX2hMUs+maKC{=lsZ6G9b{p%b)`O> zlXY~<*!l|t8j)INw9o2#iL(M!$6poRbnv zCOH7VTqfh@1N%K>HY)Vw&x5-To!InUPl5fHq-)4Qbd2hT&Fl`(0`H3*x&61*8PzdA zhkMig+tNeQZ&fw=yWs#?SjDPKzBfSMp4Ps4TZ}XOe&FW`XvXn=qJ0?{<9v%{bMIH|hbD8`K`z{TdKxt5jV_65X=@M%B{dOnlGTRE7v1clfX^H%-x817TxL6~P&#m@VLO?D?vgkuC_NNwPB+Qk>o8EeLrwuCOt34} zA)3IgqpOzy_h9tT;q7rc%7D|&$vjifajq5600w}!Vpazj<4$#kWAC!Cb$&yHK0^JCJf0sk+LI2-GU4> z8H5V~F2fT{xncl+7_KIv(%;cWCz@zvEn6G#>an5GZ?XmlG4H;eu-DBk zEn5fkus%T^CcPX&o0@SqJkwBihV% zT)1;77iO{{j|HYk9zNr_pqkyjn6?5!xI=nSsA5~MZKfQSVs(xTO3ZF_jKgpdNLgN# zEGOb~p0f4FuXjNY1J(QGAYOg~7L{4lC9INtd%*fYK3L!3B@a0gz z0gLH=%cJ=KOz~;+$kB|RsLGv*c9nv4Yt3%*!o+~!+!!|Z5UYyBAfe7v&Gd#d>l|b5zK=({HD@%A zhZ2tF=MY$e>JffW0I|QT2laU>onQ^(qU{3R<7o_Suh0v}GgC@I<@gPt!iqt3rIk*b z8+GQOEr1D>0&U3pnI!AwGXCG~c-mhb$(f0Y9eJWcKx^ajGJ)lZSrYKFa$RW6_sd@w ztQ%x-6bMKlz=*n(^K zs-tvIF7C!I?vCWGBZr=7&ZR61u$~Ojoyvi$+4_mfe+LX~fD5AVf5H$DjGc6e`xhx` zMFZJ8`8{CdT}G-CsFDG?d&dQ^!AC&Ar*~yHRH*(7WPqfixZ3b<@EOsA)qz1DWX2!; zTkb*h?M?7^N#>m{wkl9=hz(<{3LLb^lu|ENr}R&#JYX5^M>Q7Mpi zXEMipSXxBLd6-C_a(@vb5*gKSB`tANfV|(v+Zx}CF7N+M*5`!8mJ6?L@T;EE=Y=3h z$L3-1V4Byp6qVh6m-?w^1P#6>cnV(Pyupd~T$yCbB-7b`dJtBp*NL5Aac*vI;oG1S zSUy%G91vnAFUV-}C~{`Jw5?7(n@OOYc4)cpY|?WsE!S$gSwLVivpR*6^1$Sm{M?b` zeyst)yjO1$`dfC#XRlsDH55zVsHMPrImJ87k5^9mD-a6Vi)8(nQFVfSdd+PM2hYJ9o)3_`i2N0=(K>$5N^V z7Pq9!(oB9CJ(l+C*cIN(4lkM{C%v5=di1Na9PGb*;h?em>LEdb)VDU5PE7t(5-c_m z?@c|=TVWH>Ur133u&tv&A-tn=^JXB&YoAL;mzy3F?zfv@2N`Im?b+e-+W&fi`w_m5 zFMrw>*2jDuzs$IEpgc06lVPY+FC%B5J$=!|LMtg$Yr?w$@CPPZ8 zD)ZtZ!eM>fv0se7y1&=9ayd^{fO|W?NPR{LPgu^DRfz@$=cT!hic`imre*Y;j~}dS zrGDF+m@PKqUUBZ^x8g3V(nbosENvaeNgjKR$7<@DmV2r_8a;Y)31uW-gJHr&elKQR znLjZ6A3`xV#y5)Vaq^S+!CjG7H2I_K#On&pj@6BYCntlQJ?gk5o7_4FwnXM69c&{P zk0)3iP#Udm9xNE>^vX|?a^w1nmlv*J9z{8vwAh>_2|?wNqh<}3GcCP`>Xv*-%r;G` ziImjUifl`dcUr~Oa;vw`m{h3ie&uctZ<6k4=zRj=r5O&B#qe~4LHs~}pMP(PfzYv4 zMOx`oXRSK^>^&{nMmb$kQ5`7Rkg+s%MMdyN9TF?335uP1N}Gpy;8R5E78pY9!7I7U zU)3u~vU|?xywYwyIx164-Lg|sMRjH}Gl|zqd^B@QWwE}^@@@MjrLe((l z8*~Ma3?FHcArND!Ugjdszh3Q}RxIFroVYrYj!Aoblw``87LzO8HQ_SJy$+sMwBNts z2D{S%VUj`LmgD3tg=eO1EDV-4j=T5iQUqpR@h}oRZc&wadv0Z>a1j)rBAgx(-J7g; z^si1DMarBe2eJ0z?TcE1J;pl<{*Yn&3mE()Xw-7%$@E6=ZMD=p*TLoEJC3WjWTu50 zX}tq~Fn{E&IPvCSzq1zv)jG&>%j0TdBH#nYP=1-2VbR11nSb__26$F|0h(=p%q1-1 zpG3E|7U}=x9~(PE?{m>QwB9J>SgV5Z^;3(%5eD^-ifUsd(P)!smnLs>(zfW9HH}nJ z=k7}5fycU9_&TGn^EZ{Q+&;E7n!A{LYlRQvLEWM+YdLI+g@t*_baj#jn=+{Pvg>hh z1A(e`EIdsAaH#hA+C@~kpJnCJum7-u$Or98J?~(xRYWxZqnFX?g-(UP844Xu0<#xA5G8GC_)1;(|hS;}#qmzEYwZ50+u zvpl6#nCmFXe^L4BnIVMmJA|T}w7N4!PK0y`46|=6>OFa`#m?4G%*j;;~_&>QO<-#%x_B+Fur_ zI|-CHu@u^}2qDgcZr{@I4`8^ReDy-Bpz6kDk7cFHR0Q=hr#y#rP4;&*W;Rr(IZ7!y zd}+=v;1TkpmX(4Lu_0SUspk52>l5Dn4$ewBPJFD-2u3Q{Z7pI~-uS@kaD``W)mgUm z{m<}mGi#j=?Ut+^&O3j0=;xy7YR%WfzX5h%p)GH*@)Lg3m4SD`{F zjinL2cV+BnPJw-i4Ehorj_veVTGqflO>eX*PS$>Lqm#uj^`T*Yi$dT=uDIsfqQi@wyO;63-_p4>6p^^yeQvcMj zSM(804`@5<)!0QjxCC4svyKI{Gbg#)1$hdmkG4vg={(lWOgu*Odm+7Nv;UqaPScA+ z>JnP4#qTokB#$vZP;usRPi2Tl^va(8nmyZ>vfRs8T`Vu=zQZ=)^~+OEB+J}(P**w` zLW!=;aQxk#Wy`yXBD##AXf?4#jX)UajnywNxx9be_q3`(0Nd=4UyD%4I-;MDP+M-w z(4dt8b;GLsT^r>=XyHsVSXP?a&NaZkU+fsAXW@HKj>m7G9QP$tu^N;1kewwy5@EiV zO)eZvg`1=5V|3E#U=A~yhkiz=^|W%v3Iuu&Aw&Y;wwcp&RLcKRK3Zv;Z+6g zmi|9e2C}_4SG;V+2Hl!fyu1p*%mbY}Q%fI3Syvaa+tZv)kL5E`)V<-C)Y|e&gX^j0 zBwcp7ox~}3?TO`mBO9`Hf5xr{c&B)BCXY+af-s2)v znJ@)>cdUlLmt8@}VhA~|&g1(k8GLe!bDy~q=o1NEgtz}L&9OY+B_(kS7Jm@5&Sf!# z1w=51CiuH3Fu)uRy1{S%nd!@{`{gtH>M=t7KXRde1}blZROJJ1Z0^YbR=%Zn=V0)^ z;T}Z;%#kbIw}20v9T4El=fd{`q=dgd@4pRIB9ZTWGTxEtcY3sy0{K3Odq2WMWhyJ_ z8o0mDoBcP#`)A<$*9Ojl=R0m2zw7jTfUtBxd7e`I*d*u{a>)OWa&&MK5X|r-_~+rR z*Sg62}o#fE`vM@XG0p1={rZdXDY%Wfh@I*l|%+RTbp^yHH|_q#HJ{(oDeQKVfdhh<`wN zLK%oazsBwdNWCy2B`pP`+h&S9pv%DY^=dWTP@s@ed&LCWOB9>P-IKJTuv@x84DEyB zGR>qAo5H|>$SWSlb&%u`;BhJQN(pj z4Lfm)5f#j+isv51&;ATdg5$-tWJ0g7`YZB)JC;ayRLYsc<)T)1!n`>w5nG`_lsGUFYNJ<;0P z?W{6kXbWWp2-}L;l-yGg-7w;*4dC+74MO$+P~2s92tYRh6v(4w5DY)IV6_MuR!r2! z$ZIkOp{jQ!h^J4-kn};3%bc0+m`TWHL03hVH~^mLbs`Ip5~pn90>gQ}!{tWDkdT`fR!KKg&mJaq?6o>Xp91mvX|4>j(h8Atdt!7kr!| zKA4ml#taZN;Gur}v%7*G*gJ7Tz8!Sn64>Fxl3I_2sV9J7W{<%%gI;YiTsoHwEf3HP zL!&hhL|Z0cX%5x83!pB=r*}oaU=XVAx`yXI#g2VwF=LnGSp5Y3*U`b`?bA>!|6Br) zR1O*D2F?+6n~S8SEAApqEHbQ*n#i_n0Xyv72k2fKxd3!@sQ8w8pRP)+lOYRZOFx;h zRxIIwmYWu!FB@y`+`&!!Rm23KtTxh$!I2Y6Wd^UgyL%>&nYsoZJk?~;^_eE_6err1HGnWqS00V>@1ta$^_5Oeg_ zpH#=s?3v+$55XAvc4lVglq<}68UUJ(uP0Fn-t`&_3i_M<6|ja~M&Ar*dU5b5AE$$; z66Kgn8&)k4!4h;$0dVANfN>mpI*Qa&iclu%Ew9YpV4g|X(eAh&4@mn}X-;|yXps#H zY*#vHo%z7MI@u!0pHY~`<`4<*;()Mg8%QGKSIk1$>dhKxr{CoYOOf$(6g0M%fP%<0 zG}B%MZJP9K+Mw4b2olEWWo2d9w=P$M1dD*LL9oNaghg`L980whFDPIOUlXVY%ZD3bxOXtBTU+xV%K0Q7vlolq{7!{rwR~d)@{?--p$!_x z(lSZIxM(2akJ+}Rxltu#9LpVB(;vtvH_YpOYa2neY22D^J4W^kezC?f`9#`nRa&1* zq4J`a(AyUdxvke^gkL(jXMITOSYAVJuFDl%#PwaL0j>2ywX0!JMNcD$V;gj{0K`51 zdjTjP#KARWM_n+q`YClNcj{=ILRIrzAF4kxL`9T!bv;xRV~`M!%MOE920>wrY;uVg zufUtaR`|LMQDv~(h-In1DVIEk#9LuhIHnmPXX11G|utE#G_H z^AxDyKz({S>dtI$k0p{|tK!G#cH&t_*3w z=GYBYYHG*y#4w1AsQ0(u9Eggp%kVZw9nUJJ0=vqtTiwt2vB!apTIM{&nWFx=g+pOrtgy&26t0poDB+UULDvMoa4>^l zLaC4nr*Z<$B!4gDInZ!J8FT|)*jTzGt8gr70$;+rZwA!kn~*iC3^Rqe!ynr$Ff*YG z-&8yHn85}bZ|FtSBs9E&5<^iR3_O@QG>;u1hSj(Ar@W4#G+Tz^nmYQ3$WCempPv2L zl^P4gUYddNi5G+jtJps%+JS19&}a+NoJ}b18K*muM-^_rDT>fWq~ zln{RkU5Y>rk=!triyZ`^85vfZ(w9V`S{_v*`d=KZ=&?gFwQbp2Gm!d_NIuy30GF=8 z1Q8gBiNA!WpG3XPLdFu4NK6jym|EekJBbyNX!e3~X&IQ$%xVYhum}84(Iwb~Qn{ea z=#fDQ2$m&sVl*g<$**<*<0XIWFCz{FEP6o9Q7VPA+umr|E{KtSpsxUF6G!3VRK#vW zr4c%ZL36aPLED1WcbsZX*jAcxW6-5-mppW7E5nux7~4J$B}(h)=wLw*rvnWP&{}wv z#)0=DQ_tWRP>qE)CdiG}=hwrz`Ne;XeBP6xg#|%?FZl3gA94J1^9Zv0d?NHN>8QBd zMfCKUL*6)}^ikods0gIwfqXjL5uvkm|KDVZvLed$T7#eW=Rdn#L!d9n6=#19Xio^K zpZT|bkR#4UUGUU7De!k<74%(;WuzDL$MVdD#UB+NpLWsRc9&|b@_pamDoF1hJ*k~0 znDf7LG~FNs5vlj`Zxh`NQhInORX z{@~5t%@1M#Kh1shakcNY5wD|=R;q4RQ~;{m5DMz-?FE|&35eosa{50%_hzPcPdLd?zwa+$}>Lk1sN4oNqA zj#!Jher69M^8bN%XgwOaukw8mfY%E%*12IkAkyQ~2_7MO%XS?W{^!K^k+x9JZ@|^* z^B4)vufh4ZZ{-XqribqV9{>IJcfPDhlB93#k{Rzx$r~F&vqZ)IX@cn%;(;Kt*cuv>(&##to@&B z8S;M)D)x_S;(!Qjj}|zX^q{t7;b_c{>(!p~{|da}Y{Q>=@hq8bBv1I0YF{cqR8sSt zD+>4{+f8_KA8NjL&7IQ#SpD}DeP;@Gj=rERqPC|T`33$vktISPzI)15td0kb;lLC> z!C|fh0}2l22XuZi7l`ubajQ>1I9p?}ux8=Y@NjZI!orK~)9QI^a`$O-Ehsnov?x9; zinm9_r$sSO8U6oseRx@qkrxpN=A;Va@yAhVv@ph{MN=%V=rRLH0TV(xe3gC7$^=qU zd@paLPX!i_Tb9PxU;p0Zv@e-`@it98PbEi}^`ND`B;ItitI}m>`JbjqV`ODmGqK)B z!}_4;4LQu35LK(-Qqr5;~5mRt=y&EY2?P0Cw)ro*?R zO?u<7i@W;ntPCrA{kWA(|EUWl2f9^CZ*^Z-IdK7={8lA2&{ddc=ue_dD<$exKaJ{h z6lmTlKiO2|T`;~@FA{`8Po6ffQYbRjX?3SnTy|h~xbg@RzniArSXNraZ(@|$xPyQ*hAOnq2euc_J&z>htoNF;mPEWOj^r!3N)UopKVCBM+!ht*$` zE?y}0;>n&8Ln#>bm9i+48f(I!w!C;84|l~%;G7lm-c-jblkeJ} zKGQmp5R*_Oo~rL+QvE2q(|=%Q(-kkT7(>0USaRZ7pw0j%_e>BQ_!j}5Nv7e#jeO8v@0wl3hYWm-;8VtT1Y+rMwI~ zU+coWwh7cWfXO5w-w4eMDOTF+9bUg`-*j-(jPr8hZ38YxkvmlhqsLbsMMDJY1%>ro z@N@Nvr$rLAa#~NE5Ho5yV!gF@hMlxzS{??&Qw$xpK&vAl@I}Hk)1>oti>ZQ5!GS*v z>~<6xX03&XSZ>4fLNb$X?Xf&A-Lf?Q!D!+iI$%yv8j`|y)U0tw!SXnaQFUikSs5#u zbDW^prd;3-tK`$6MbHnLd9LDAE$VxM^Ncjz>3utE^Q*3Ra#)LN7q%q}p zcIeTwvvS>!zZp+{*nV-4xTmbh>I8^_@ZGSsPdo1*2J3|%Z@wz3@;eF3CH>T-qDRVKIb~5GA!7kuD$zg`#b+i zz&%7V-&)?F;NtvcQ`~XQagv27d}VuL@MzQ&{BFeOQL?mA*COauK=cP;D)Wl)l=w93 zmbPxbN{*hlguMGtH=Rewn(h+c@{O)e@lKg}PBH0CP|b1eUP zk;4~pU_P0W=qgHf6%E{J^Wy25JRXK{_?tnZ$*b%8BwZ!GT_TobWjjPs5H2dD0}+`# zMPd+nS|z8Mug)186$weY67f-joJXZwz8J`6`v-Lm6AgUNFW0iR8|oK7*OmRvxc5eX z{fT~NH2@RI^`tb69v4p?gA}F_=yU8hNamQEi)Oc0g z`E90F+K+Fh?KySG5W24J@V_Rads}0>GR%hc+IqyC*(1*nh}wj>!K#BK(fuB zcoZTIf28{9(^YN!2gOewW;BulrqgI%XJzmf`Dv~L1Ejh=u0`;tIREMvu>131&0?ds z-(6U9{^!V NfU>$$@()J8{U5L$6X*Z{ literal 0 HcmV?d00001 diff --git a/documentation/gis.rst b/documentation/gis.rst index e57d6886d..cdde15703 100644 --- a/documentation/gis.rst +++ b/documentation/gis.rst @@ -31,6 +31,12 @@ >>> hydrant_data = gpd.read_file(examples_dir+'/data/Net1_hydrant_data.geojson') >>> valve_data = gpd.read_file(examples_dir+'/data/Net1_valve_data.geojson') +.. doctest:: + :hide: + :skipif: gpd is None or rasterio is None + + >>> elevation_data_path = examples_dir+'/data/Net1_elevation_data.tif' + .. _geospatial: Geospatial capabilities @@ -47,7 +53,8 @@ The following section describes capabilities in WTNR that use GeoPandas GeoDataF .. note:: Functions that use GeoDataFrames require the Python package **geopandas** :cite:p:`jvfm21` - and **rtree** :cite:p:`rtree`. Both are optional dependencies of WNTR. + and **rtree** :cite:p:`rtree`, and functions that use raster files require the + Python package **rasterio**. All three are optional dependencies of WNTR. Note that **shapely** is installed with geopandas. The following examples use a water network generated from Net1.inp. @@ -822,3 +829,101 @@ the census tracts (polygons) is different than the junction and pipe attributes. :alt: Intersection of junctions and pipes with mean income demographic data in EPANET example Net1 Net1 with mean income demographic data intersected with junctions and pipes. + +Sample raster at points geometries +-------------------------------------- + +The :class:`~wntr.gis.sample_raster` function can be used to sample a raster file at point geometries, +such as the nodes of a water network. A common use case for this function is to assign elevation to the +nodes of a water network model, however other geospatial information such as climate or hazard data could be sampled +using this function. + +The network file, Net1.inp, in EPSG:4326 CRS is used in the example below. +The raster data in the GeoTIFF format is also in EPSG:4326 CRS. +See :ref:`crs` for more information. + +.. doctest:: + :skipif: gpd is None + + >>> wn = wntr.network.WaterNetworkModel('networks/Net1.inp') # doctest: +SKIP + >>> wn_gis = wntr.network.to_gis(wn, crs='EPSG:4326') + +Sample elevations at junctions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Elevation is an essential attribute for accurate simulation of pressure in a water network and is +commonly provided in GeoTIFF (.tif) files. The following example shows how such files can be sampled +and assigned to the junctions and tanks of a network. Note that elevation data generally needs +to be adjusted to account for buried pipes. + +.. doctest:: + :skipif: gpd is None or rasterio is None + + >>> elevation_data_path = 'data/Net1_elevation_data.tif' # doctest: +SKIP + >>> junctions = wn_gis.junctions + >>> junction_elevations = wntr.gis.sample_raster(junctions, elevation_data_path) + >>> print(junction_elevations) + name + 10 1400.0 + 11 2100.0 + 12 3500.0 + 13 4900.0 + 21 1200.0 + 22 2000.0 + 23 2800.0 + 31 300.0 + 32 500.0 + dtype: float64 + +.. doctest:: + :skipif: gpd is None or rasterio is None + + >>> tanks = wn_gis.tanks + >>> tank_elevations = wntr.gis.sample_raster(tanks, elevation_data_path) + >>> print(tank_elevations) + name + 2 4500.0 + dtype: float64 + +To use these elevations for hydraulic simulations, +they need to be added to the water network object. + +.. doctest:: + :skipif: gpd is None or rasterio is None + + >>> for junction_name in wn.junction_name_list: + ... junction = wn.get_node(junction_name) + ... junction.elevation = junction_elevations[junction_name] + +.. doctest:: + :skipif: gpd is None or rasterio is None + + >>> for tank_name in wn.tank_name_list: + ... tank = wn.get_node(tank_name) + ... tank.elevation = tank_elevations[tank_name] + +The sampled elevations can be plotted as follows. The +resulting :numref:`fig-sample-elevations` illustrates Net1 with the elevations +sampled from the raster file. + +.. doctest:: + :skipif: gpd is None or rasterio is None + + >>> ax = wntr.graphics.plot_network(wn, node_attribute="elevation", link_width=1.5, + ... node_size=40, node_colorbar_label='Raster Elevation') + +.. doctest:: + :skipif: gpd is None or rasterio is None + :hide: + + >>> bounds = ax.axis('equal') + >>> plt.tight_layout() + >>> plt.savefig('sample_elevations.png', dpi=300) + >>> plt.close() + +.. _fig-sample-elevations: +.. figure:: figures/sample_elevations.png + :width: 640 + :alt: Net1 with elevations sampled from raster. + + Net1 with elevations sampled from raster. diff --git a/documentation/installation.rst b/documentation/installation.rst index d32013c43..4d6b8bffb 100644 --- a/documentation/installation.rst +++ b/documentation/installation.rst @@ -272,6 +272,8 @@ The following Python packages are optional: https://pypi.org/project/utm/ * geopandas :cite:p:`jvfm21`: used to work with geospatial data, https://geopandas.org/ +* rasterio :cite:p:`rasterio`: used to work with raster data, + https://rasterio.readthedocs.io/ * rtree :cite:p:`rtree`: used for overlay operations in geopandas, https://rtree.readthedocs.io/ * openpyxl :cite:p:`gacl18`: used to read/write to Microsoft® Excel® spreadsheets, diff --git a/documentation/references.bib b/documentation/references.bib index 40f09982c..50df2e953 100644 --- a/documentation/references.bib +++ b/documentation/references.bib @@ -128,6 +128,14 @@ @misc{jvfm21 year = "2021" } +@misc{rasterio, + author = {Sean Gillies and others}, + organization = {Mapbox}, + title = {{Rasterio: geospatial raster I/O for {Python} programmers}}, + year = {2013--}, + url = {"https://github.com/rasterio/rasterio"} +} + @article{jcmg11, author = "Joyner, David and \v{C}ert\'{i}k, Ond\v{r}ej and Meurer, Aaron and Granger, Brian E.", address = "New York, NY, USA", diff --git a/examples/data/Net1_elevation_data.tif b/examples/data/Net1_elevation_data.tif new file mode 100644 index 0000000000000000000000000000000000000000..6601807069242bcec505890ff7b93a2ee0f562a1 GIT binary patch literal 33444 zcmZ{s0c@P-dEZ~#ntMq(_v)!G;j!sgM|BB7@D48P9t^>KqiGqFYQZ(-T+Y+fTfk{)?|Ne$&Q;t?s(D;0c;5omBMd^gx@tfeg7?87TsCFh zgC=C0ry1k<`~8IA1IEJm^E=y$%wPY{^M9WA{q&KK><^U?LN$b2ZVR{Ue)p5Rf8E0O zKezk+t$ZJMzrT&|hyUE|-}C(s-m)kBt=-QycK`aJ@A>>LzTdn1`48`YKm2Dd@BTIR zJ^#$l^ZjkRpWpGhkucQl-M>!HkA&NIbueDT@t5Av4qR`r(3*Pq(`yW#fto}c;P`zpPshX2m@ zzo#KAeCov6A6b5E`0K^zr=~yj3;*!vllsJ|-N)UpEOdY_|6%cbpjYUVXpTPk_2TnI^lfzV8^z~c^alFy9~GbX(Dz>`o@daP(b4}@ zeEt}^j^6uH@p+2Aj!yjN;`3+FE%eZ8@%cHl|GUNW1YDNja^v3@pZDPM3ApU-7N5U@ z&cS6rT&|+ez$L=v8oCUZhv8Cz%Qw;EaQOgS`siu6Ou*$WbO|mWh06`}7vRxD-v@_J zpnr_+g~O-NSI{4Y#~k{n=ze&74!w&06g-|m{{oHh_zHRr{aJV{qyHCu7%nMXa(olv z(t%3~ml|Ahv<;UCmjT*=O9L(e9zC=PmnvK`v;~(0E`2nGOARhL+J;Lk`KAMxad@QY z23mv10-B?Bc(lQWspZMbaVn*<)~=p0=3!KII$g3CC5xq_x}IRKXpbRI4>xLidS;4%T196bY< zI$SnOzOmnI!DkUJ2^=0r*U@|6F^678_rc>N+DGq&$0_tRbQ~T}pjXiS@JP{bqX*z| z8r?wegU3ADhRY_tsl(+ATylK#EL;x50NJMt=-0C(!?n-VKjWp&9yj;P6TGzn~*<_(k*^Xa_D^_@)MzIk;r_ zW)D18&?9gepo`riUoxr(N6*$0<3^f+9q@VJ0Z!{r9rfy?2NZ|pbbH0PT+xQxLi zLucT!2M!OTE9g7u5qLa@4$wpJSVZ4K>+pCQy^7ukj}(0k-3O1y(KYmLcpOI?a9M=Q z2EG}COA{{3_~ry$?t+K)?2&F9VhK08aVX(AxJ0q~23l9GRwXVGn4(7E#;MqC!AHm}}^l#7z zkKaMxM1K+<^XQ+W&Z7Y?_8I4y4BrH}M7Tt_xCcn#lXT+{N65NyD4_|L0FMr84#?r- zJ|KZhhBo06;Gu8qH|8|w8}|Yke420xaOt2CKIVZ4k2dOlz&R&{PajR-(Lyt{3Wp|Y z-wAMNpdEC8nhV-+$RAOY5iWCZ>6Pl4`OJB9 z5-yu)3YRf>tfEtJ8K7t2G7gtZCEwU@%xTUy?giYZRp7FWMz~x>=isszEk0rDQmo>BrpB=Pa@{Rq*oaTJv zULb*sd*}vSw$Kzlqj2e=lW@t=Ik;5dF@<)}ZL|)LdDOl$28UVHc_+Z(2V=P+{^(F!^M5TX}EYU@FH9e!{v3f zg>P=hFP;zFk6+f%WAM3-J`ET1z)SEs0vGoK*0np}VqW_ZI5>y?5zp?|J^~ld4gLqu z`{D9ysJY<#;o$kf@A2%u;HTkn8cpCb36~aJmf(`%n=SZM@k^b0R>Los&}sa#ji%i= z#1r`S;05!@EO48 zX}FBR<3%)fi(mgWbOtWpMHk?5H(bu658;^lQSb~f5&Nv)esPkxmgSpM~oee%u;4q6i_tfBW1TM31 znTJaU-&};x2EMtDZ>sp^F8s2L-Vc`y{TN)XqL0F72YnLX?8Pr#w1!_^MvuVdb#xX! z-$m!)awl9m=zZ|Hh<+3<8|WwDaveR5Z>sp^IrJ`g{4;bJ{VVi-xcoYrq2C9GkD*^e zuc9A>%cE%ir5J`^JMb|#kHf`1z%*Q1aPdAshkCXOA9KKU__z-k>BiyzOjzp1p@e-U z-`H==Y0fw91>C2(hxYu&9Bm#Lhsy=j{lGeWrs3jwK?^Q7Q1{a2Y4-#N;IfL|0|#@# zi#%UPy-(o2;1&3ot3M78&k^3>c^_Qba9M`SD!y4S)id*%^QQYV&jUQ)^uESRaEbWk zRdlKwhgiZJ-8hsm3!iJ~DSUG~ercn3Oo{=)-Wau00AD^V$>eaSr_= zT->ick8eCT_$L1GJoGYryyx*5c)W!EA2h<_^XRMS55Qpx{R%3d2$u$2ns8~sCB-+M zpLOs}58q_evp&Ab@y#4u0z4MbYBvsX-1Lq8#+>GS<6a=b#XYp=H=cVo;o|vP3obio z3ZDvGJZBq)O9vf?i+#uQwhTVbJAL>}!^J(@1RUm2&vybGJb&|E!8)H;OMLXtV6pf) z-GEO8E_>mUz~umarVf|GaB0BhC|qXY@*rF?^uusjM?c;z{&zRfWB6qgeH7me(8u9( z13iszDsXuQ-3yOr(F88fqX*#g4YUrIm(j!Uc@=HI<=f~{xV(wZ!sl)DLAd-TnxVJB z;lt<`Q2FQ|>stpdD{$$R>Y4e>dDH!w=K=0%yst5U&klT2xa=wUW&|E>)cMA}fcrG} z(4OD8w>bcp9y$RR>)JuMnAhCjIENmAi~F@HxOi@G6fW*@X5ivIkK=GLx1E5?4(dK9 zg^T;PJtaQ+$NFafX~1O)F4J&1if`O&&fuF_d~>{%15V%@_nW8S65;YK)cXOiz-baL zucAlr%Vo5IUtUAa1K-9s({OP=a1=gopq>}J38xu+S|cwQ6W65ye4 z>^J7L0FM-PpO(VK^P3bd0Ujyp{WRyC0FM+6@JP`B7yC|tM~XV{q;LuFNYMa~6bAT_Mxxrm< z@jTT1wDag4aPj`y`{81)z8x;}sOO{d(LdHV`;Ym+{ARvx!etsR4SX|+Zz8^_118`YF?=11rKh`(Uk zSB84e#=7RcPV<`Q1m$w&WK-|Ro; z1M{2tI;S5wpSnM2!R08vnZz#>_$9$FBlu+p-NZM2d~*@stl*nPe3RmvdKBNhhEC#_m(dCQ@&cN`=PWvcUrwVt=%?Ybi9UjF`sj!8%|-P7 z5+D6zeY5|V56o}oYv&*5Q}+k%cM`Z%;4(lr;FG~;1wL)K%)zAzmr1zj8~cqp&H2W? zfcrG}(4OCT?&*Dj1Rjg1_d2}~J%Eq*)i>aC1TGmm0hblj`|EAE?19G|x{WsB;(hkX z5+D6zeY5|V56o}oYv&*5Q}+k%cUJLF|94^-eqDgi5?oHfr3shAa2bco2wXgWz5<_B z_?*KxDZZJ(H)N4i9R4;BGvK+|jYB-4-;F~F5kB^v1RibF zd8Y-J9Cgo@!KVfnbG3VF=N|d!AM2a_$9!OZGhaLZIG^U^1ot}`d7(y*>!3BbjKF0J zy#ya~z$v)6517C=o(pWl$9>H*zDe=TQMmZNfcI<61MB!_3NG#kYVddt9qGm)j&QOY zhZ45BaVWw21DD`*KU~ZO-oH5o7vC>%U+@Zi+~;)RBp>}_eY5|V56o}oYv&*5Q}+k% zcig|ZA9a82`9j_<{<~J;(}v3|Tq3?1h06e4EBVHLV@`9vaWCLL%{{d9XOp^CFV(dZ z@G-9)fs1qKzHagRZlgC)&ke4?C&I^h)ZFHM+yR`9!o_>J^3gxmH~WwI!2D*ucK&fb z4des&JMQ1ykGjA1e8C*ohRY0G>Tnr>OO9TE&jMVI!lj09_TZNd^gR70#Wz#PHaW%$g&=U%w9(Ra`p_*_Qo@aUo=@OT`}(T~FC0{Q`XETGp)eDsg?&HiIPFu$3v zoqwE9-5N-6eDsg?&HiIPFu$3voqwE9 zYvcp>JAwS>e$@T7=L;SD)8(qaW^YHOJZ4BRRp%?It z=W0jr%_v-aKh3&!s+8CE;v47ChvDLW?H+h|-sbyjQ{6cH&k5h`#-W4+9#5j)U;8+| zk&phdzS)1w2j(~Pweyejsrv)>JMQ1ykGjA1e8KadQMh=1>iKqr%Pd^X0Xy()P3|6{bm#{S5foy41D(Do5xY_ zN4uXs2bXtIbM?2n#eb)K^pEw;{$oBcznQO{f1FRRWJU^|# z#rp^AJO{X>_$Gx*giD~F=^OivInDVdgOB?(_t1fQmX+$7bB=i}!#~cU_MJe#$x!E= z1THCRZj12fpaCB8$#}ND*?)3;WPUSWJO4E4LmBFRCnd+_sQYVkTySrIW@v(staCZf z&O^=LF8=>$13u$$*+e_=ISQXq{Bi|t({GNzrGjrRp((yO2$vo70(q^4Zw|m^3th%P zO?)#3mko3g|BS=sNpurERN|w5tZ(0hkNLp-X1;d*aXxi_;C{#boBL7s*PbtU{*!*S z`1yK%>iM?!4|d?>{T=UDWoL?iuLhS6zF8>w#(rZ?bG~se;6BYg^dwvYT)h9}Jsa<* z!rHp ze!#ru`KNQ}9(?0|&HDr9f-CTG9z6&j_XWORARql>eY5|V56o|qaOt8qQ1=J!cig|Z zA9a6iPVoF^9ZvVc#X9Hx1Lsxm@0jPj|7KoDe=COJR}P<8@{Rq*oaTJvULaD>+(Wy! zX}~4m8}HdP@sIl(`%VkrB=B+GN%2n&E;-t!zQqzB{bPNz|CkTVZ{};~pL#bAe+;4B zjYGM)--+}c_oMEwJzua7?cl2hTs+?n@NpjU{!Rt|G)q2m?~^}W{66Y%Y2%wM_$2sd z9i78B`|wL2Jq4d}d~*d&;d1~!8|XZIYV@0{=mLBu;FF_g;8TaoW{HpfvA)@V66%`y z&3x_r<9zD=U|+ZRd7MQ1sQYWr7d-#*e9rSz&$su(Cq=!#X{$?=ZJ*!N*)JAN^x}v;UY6 z%x~sv=O5=&_XqBGR>=z^@MxkJ&jmaykDh{yr1aYYCke}cfJtA z@M{Y`HGDINZ!&zd2fwVKNAS%6UBowa_*_L(`0T?sYv^%&Q-#X~bQ&%<&<_4N3>W$6 zAM2a_$9!OZGhaLZIG?&dSiw*3-`tP7zaGHpP&W>LOkuGbhZ5fE#-ZH&b0oZfa20Oe z-$~)*{WtF?dVh2cKY9P!`|$xT`o?}^PIJC-FA(8k9!TMmz{NSI36}sD`%Z*Qjymro zaLLdnT-+DPNB>yg>_6s%jQVE2cK+$mha&kcMC9HG^XL%F$6 zsKTR(>Z4%aLG4Eao+Dh$bve(?h2qy=f=?gcRLlYRW)9!<$k8X@V_iEBAM@HKoSZ|) z;NyPH`<ZIZI|die zx4nO0o;%wuevX@{_urgby+7*x-bdl$e$Mw7tasT{#jiJjPXd<){U(J^4?a13Dsbtb zb-47=S-1pz zZ)$M4giga{8%^P}7d|WKLAdnMnGzrUV|}y#m=DZv=4RWJU{h(+xrLJ54i(A)2R2~wt0R4E-Csd>iu}%Us!?9HS{1{(kEjWes$oJ!6%1L zfQ$2I#5Ya&nAaL`vG26t65!&zlfcE?7E65ekM+&|V?O9m-^|y}KhCG_58Uqr_?YAD zL!K{K=RBWlmgXVzoIYyc8}C=ykG!AQf{VE>bmQ>H8?1NjZX7;_jcy#u&2tF%4ECk` z^Tprafln1K4WAaiaX-+3&no^gPhW?Vxxo7aOK@^OEg$`3eY5|V56o}oYv-RT z`QRMte#iZrIl=w4=L??yw5W5QpPIM5e_&ql{?00S!TWFKZtss?hnxEabA#_MjC6~i z_s7vC^bOQ~cZ5p=E=_#Xf=>#cHg&B7pB{WN{L_a|4xc%=1bpLpTNN(y(LdHV`;Yk` zQs2zi&OgqlP5k42rv)GPqbWJt^MyA3$n&`lej105KB}`%$WZ%HA3sgQ$6Pl7k2%!& zDBz#NsC$NWo~!VfDsk$}$1wa_flm)UtEIZO1|Rd<27KJ#Y{AF<+75hD`0Oe1(LdHV z`;Ym+{ARv({&7A%Kt6E4lY{1F;J2||(|F#7u?~l5l zI{~L1)b|%sIGOMEz{B_F^lM91h_cw1h@q2S;;@v zH~WwI!2Fidhn#;>`cS|>DeC^sJ{0gziU#0WK*T@K1^c zxY(BhTp}u`jA#86;bq-xey;G#G<+KHnS@V-PaQrrd@~N81mBFB*Wgp3uE|ILSl{eF z<^%Ja`P%uX2_N?d?swe3xgT|Z?fHV|KM_8jpVr~y{ev3*@&3*@e7ygbz{mTe-tS$6 zR~0_KzfggX@Bhfj_viG}gK)9#y~MNqDLY+o>A}ak)`pLHtpy+F&_*d2)Zk-otCsv@ zeY5|V56o}oYv&*5(~N%Pe#e~Pe$+bWeSs-v;N<;-75b9*ciQ;L`)_mb@&0HN zKY9Op5@;oiLtiq=apIP`s_>AJ40lHS= zqkpV#_8;?s`OSRo{NsFDhmZRm_ixs@BXC(oJ^wM!d44(zC+{CbIC+1^dguK&=U4BK zuEEKCCnw)u&`-YqW8L%pIs23EpPD1w-dEgSsD$_-Cc$qdEBGsQsvipGM(huG`|d0Vn691ph3c za+-#ZesZ5`-Shn{`;+fCnj_q=n)hb$lexT3pX#9QTlH1^*M(nZ;ZuXl4w}Jdp~Of3 zSl{eF<^%Ja`P%u%`PBV^Il=v#`%(AT5g&Q}V_xw5v<4UNAMEh#{T*|I`-KJgcz;w* z-oMsQzQ17I^Zg(Dlkd-&Bd$>Y%zN%1%;o<12j+!!V?9T6h0|@B3uGo zgsQJx&?fjF`huk07hpclMK1#{k`p9`G!$9;UbohrE0;M0N6LWz(5vA)@V z%m?N-^R@Gj^Qrp-_dD+2a`K$}YtI)v|H-%Ed^;F7~9miXu& z>zn*Rz!It!-&pCz;gr!}+*pB=PK-5W1)n#M=^X`xj2R^eoS z+A8IUQ8=0RCgGH$b8xD_XA140+h`p=^CezmaB{z=zXF_&pw`1po@?;2e|i3tK3@F1 zJMhWiljEZRpTzzHpJpk?HQ=L=_}`3b6+2BO*q-Ftb_7$Z?ps_{WT7!7HT~V@N&Mge{Jx20w42+bus>0!KDG8 zCVX1(N#WDRM;$oz;FQ6s52u_wHwUKxp9QpvkCsq5t(5$<3McE{8h)}rZQ!dWytdF5 zymrtOUsd4b{4|Q6I_NlcP+oO7WvKq@!)>}$4=3Q{zTEy5;B*)@Z>;lNg^&5r^TKeV z;1b{y;1l2z;A0;O@X<#BJ^?=VqX3@(pR~luz7*h-qSn0>PWGo1PQku}n)mcofKQ4# zKc#T8U&$-jr%L{^Zv}XzXn<3S1~{c?fR{cC@QP@PW~kgEoa|@WV+EJquNK;YQyWg! zxfYz7aB7tDLakKqswGZ?lAp|X*1c)`WPh53Q-;<{^=}+b<_Y%?&QBHmBro%y`w0Ej zgx3aYJ#;@|9<+be@Yf1zU0i}&6<+2;xw*fvKF*f>=KN(|Og~fbXu>DLCn@_;h_nk7yNob*#(s(YPMf2!f94r(6PR|B{;;nYWCiB|%zcFA8Y{ACWb9%k^X;V=7F z)jEipH=5MNekmXJ@LK{e`QAjwj;QJaQNWM&V^W+~V2&tNkm1*8*zZn1+|! z%!h6G?V$Q?4SqFvIe%FvXW{0&T!)|i&HaY_!lw%^>Aytb6yX!#W4)78Ecq!b)xB8i zPv$&xg!7VpD#KrCiI@EiI_ zznq`+R|>BhysU?9`d5UP{i}n|#!K8bc&@?Ae5l{*@M@!5sCCl$%YD;iiQf*-4Y;}A z2=LN>J)SH24RwEL{S3bhk4FoYQwAsf6yapu%kh={DJl7?2`BTOc_PBger2vm@K;u< zhvvXesek3vMY$#RFSJ?89riQ#7Up^DWV_@$`K55|qY1oPXogntTN8CY3-T&)9Pk|B KXkN)bAO9azG^ObP literal 0 HcmV?d00001 diff --git a/requirements.txt b/requirements.txt index da4347a42..45582cf7d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -13,6 +13,7 @@ utm openpyxl geopandas<1.0 fiona<1.10 +rasterio rtree # Documentation diff --git a/wntr/gis/__init__.py b/wntr/gis/__init__.py index babbb2afb..ead669666 100644 --- a/wntr/gis/__init__.py +++ b/wntr/gis/__init__.py @@ -3,5 +3,5 @@ and GIS formatted data and geospatial functions to snap data and find intersections. """ from wntr.gis.network import WaterNetworkGIS -from wntr.gis.geospatial import snap, intersect +from wntr.gis.geospatial import snap, intersect, sample_raster diff --git a/wntr/gis/geospatial.py b/wntr/gis/geospatial.py index 332b7af37..8b2989bc2 100644 --- a/wntr/gis/geospatial.py +++ b/wntr/gis/geospatial.py @@ -6,6 +6,7 @@ import pandas as pd import numpy as np + try: from shapely.geometry import MultiPoint, LineString, Point, shape has_shapely = True @@ -18,6 +19,13 @@ except ModuleNotFoundError: gpd = None has_geopandas = False + +try: + import rasterio + has_rasterio = True +except ModuleNotFoundError: + rasterio = None + has_rasterio = False def snap(A, B, tolerance): @@ -57,9 +65,9 @@ def snap(A, B, tolerance): if not has_shapely or not has_geopandas: raise ModuleNotFoundError('shapley and geopandas are required') - isinstance(A, gpd.GeoDataFrame) + assert isinstance(A, gpd.GeoDataFrame) assert(A['geometry'].geom_type).isin(['Point']).all() - isinstance(B, gpd.GeoDataFrame) + assert isinstance(B, gpd.GeoDataFrame) assert (B['geometry'].geom_type).isin(['Point', 'LineString', 'MultiLineString']).all() assert A.crs == B.crs @@ -197,18 +205,18 @@ def intersect(A, B, B_value=None, include_background=False, background_value=0): if not has_shapely or not has_geopandas: raise ModuleNotFoundError('shapley and geopandas are required') - isinstance(A, gpd.GeoDataFrame) + assert isinstance(A, gpd.GeoDataFrame) assert (A['geometry'].geom_type).isin(['Point', 'LineString', 'MultiLineString', 'Polygon', 'MultiPolygon']).all() - isinstance(B, gpd.GeoDataFrame) + assert isinstance(B, gpd.GeoDataFrame) assert (B['geometry'].geom_type).isin(['Point', 'LineString', 'MultiLineString', 'Polygon', 'MultiPolygon']).all() if isinstance(B_value, str): assert B_value in B.columns - isinstance(include_background, bool) - isinstance(background_value, (int, float)) + assert isinstance(include_background, bool) + assert isinstance(background_value, (int, float)) assert A.crs == B.crs, "A and B must have the same crs." if include_background: @@ -280,4 +288,50 @@ def intersect(A, B, B_value=None, include_background=False, background_value=0): stats.index.name = None - return stats \ No newline at end of file + return stats + + +def sample_raster(A, filepath, bands=1): + """Sample a raster (e.g., GeoTIFF file) using Points in GeoDataFrame A. + + This function can take either a filepath to a raster or a virtual raster + (VRT), which combines multiple raster tiles into a single object. The + function opens the raster(s) and samples it at the coordinates of the point + geometries in A. This function assigns nan to values that match the + raster's `nodata` attribute. These sampled values are returned as a Series + which has an index matching A. + + Parameters + ---------- + A : GeoDataFrame + GeoDataFrame containing Point geometries + filepath : str + Path to raster or alternatively a virtual raster (VRT) + bands : int or list[int] (optional, default = 1) + Index or indices of the bands to sample (using 1-based indexing) + + Returns + ------- + Series + Pandas Series containing the sampled values for each geometry in gdf + """ + # further functionality could include other geometries (Line, Polygon), + # and use of multiprocessing to speed up querying. + if not has_rasterio: + raise ModuleNotFoundError('rasterio is required') + + assert (A['geometry'].geom_type == "Point").all() + assert isinstance(filepath, str) + assert isinstance(bands, (int, list)) + + with rasterio.open(filepath) as raster: + xys = zip(A.geometry.x, A.geometry.y) + + values = np.array( + tuple(raster.sample(xys, bands)), dtype=float # force to float to allow for conversion of nodata to nan + ).squeeze() + + values[values == raster.nodata] = np.nan + values = pd.Series(values, index=A.index) + + return values diff --git a/wntr/tests/test_gis.py b/wntr/tests/test_gis.py index 31513f265..0872e179f 100644 --- a/wntr/tests/test_gis.py +++ b/wntr/tests/test_gis.py @@ -23,6 +23,13 @@ gpd = None has_geopandas = False +try: + import rasterio + has_rasterio = True +except ModuleNotFoundError: + rasterio = None + has_rasterio = False + testdir = dirname(abspath(str(__file__))) datadir = join(testdir, "networks_for_testing") ex_datadir = join(testdir, "..", "..", "examples", "networks") @@ -69,7 +76,7 @@ def setUpClass(self): df = pd.DataFrame(point_data) self.points = gpd.GeoDataFrame(df, crs=None) - + @classmethod def tearDownClass(self): pass @@ -311,5 +318,48 @@ def test_snap_points_to_lines(self): assert_frame_equal(pd.DataFrame(snapped_points), expected, check_dtype=False) +@unittest.skipIf(not has_rasterio, + "Cannot test raster capabilities: rasterio is missing") +class TestRaster(unittest.TestCase): + @classmethod + def setUpClass(self): + # use net1 junctions as example points + inp_file = join(ex_datadir, "Net1.inp") + wn = wntr.network.WaterNetworkModel(inp_file) + wn_gis = wn.to_gis(crs="EPSG:4326") + points = pd.concat((wn_gis.junctions, wn_gis.tanks)) + self.points = points + + min_lon, min_lat, max_lon, max_lat = self.points.total_bounds + + resolution = 1.0 + + # adjust to include boundary + max_lon += resolution + min_lat -= resolution + + lon_values = np.arange(min_lon, max_lon, resolution) + lat_values = np.arange(max_lat, min_lat, -resolution) # Decreasing order for latitudes + raster_data = np.outer(lat_values,lon_values) # value is product of coordinate + + transform = rasterio.transform.from_origin(min_lon, max_lat, resolution, resolution) + with rasterio.open( + "test_raster.tif", "w", driver="GTiff", height=raster_data.shape[0], width=raster_data.shape[1], + count=1, dtype=raster_data.dtype, crs="EPSG:4326", transform=transform) as file: + file.write(raster_data, 1) + + @classmethod + def tearDownClass(self): + pass + + def test_sample_raster(self): + raster_values = wntr.gis.sample_raster(self.points, "test_raster.tif") + assert (raster_values.index == self.points.index).all() + + # values should be product of coordinates + expected_values = self.points.apply(lambda row: row.geometry.x * row.geometry.y, axis=1) + assert np.isclose(raster_values.values, expected_values, atol=1e-5).all() + + if __name__ == "__main__": unittest.main() \ No newline at end of file From 26b433a110c10b877e162f483a1b130246a95466 Mon Sep 17 00:00:00 2001 From: Katherine Klise Date: Tue, 12 Nov 2024 09:20:53 -0800 Subject: [PATCH 4/5] Bug fix in valid GIS names used to create water network models (#452) * bug fix, missing commas in _base_attributes * Updates to remove the use of node_type and link_type in WaterNetworkGIS GeoDataFrames * updated/added tests * removed user tests, geodataframes no longer store node_type * Added "name" back into the list of valid gis columns, updated tests * updated documentation * minor updates --------- Co-authored-by: kbonney --- documentation/gis.rst | 42 ++++++++++++------------ documentation/model_io.rst | 39 ++++++++++++---------- wntr/gis/network.py | 67 ++++++++++++++++++++++++++++++-------- wntr/network/elements.py | 4 +-- wntr/network/io.py | 5 +-- wntr/tests/test_gis.py | 30 +++++++++-------- wntr/tests/test_network.py | 26 ++++++++++++++- 7 files changed, 141 insertions(+), 72 deletions(-) diff --git a/documentation/gis.rst b/documentation/gis.rst index cdde15703..6a4428e84 100644 --- a/documentation/gis.rst +++ b/documentation/gis.rst @@ -119,13 +119,13 @@ For example, the junctions GeoDataFrame contains the following information: :skipif: gpd is None >>> print(wn_gis.junctions.head()) - node_type elevation initial_quality geometry - name - 10 Junction 216.408 5.000e-04 POINT (20.00000 70.00000) - 11 Junction 216.408 5.000e-04 POINT (30.00000 70.00000) - 12 Junction 213.360 5.000e-04 POINT (50.00000 70.00000) - 13 Junction 211.836 5.000e-04 POINT (70.00000 70.00000) - 21 Junction 213.360 5.000e-04 POINT (30.00000 40.00000) + elevation initial_quality geometry + name + 10 216.408 5.000e-04 POINT (20.00000 70.00000) + 11 216.408 5.000e-04 POINT (30.00000 70.00000) + 12 213.360 5.000e-04 POINT (50.00000 70.00000) + 13 211.836 5.000e-04 POINT (70.00000 70.00000) + 21 213.360 5.000e-04 POINT (30.00000 40.00000) Each GeoDataFrame contains attributes and geometry: @@ -341,23 +341,23 @@ and then translates the GeoDataFrames coordinates to EPSG:3857. >>> wn_gis = wntr.network.to_gis(wn, crs='EPSG:4326') >>> print(wn_gis.junctions.head()) - node_type elevation initial_quality geometry - name - 10 Junction 216.408 5.000e-04 POINT (20.00000 70.00000) - 11 Junction 216.408 5.000e-04 POINT (30.00000 70.00000) - 12 Junction 213.360 5.000e-04 POINT (50.00000 70.00000) - 13 Junction 211.836 5.000e-04 POINT (70.00000 70.00000) - 21 Junction 213.360 5.000e-04 POINT (30.00000 40.00000) + elevation initial_quality geometry + name + 10 216.408 5.000e-04 POINT (20.00000 70.00000) + 11 216.408 5.000e-04 POINT (30.00000 70.00000) + 12 213.360 5.000e-04 POINT (50.00000 70.00000) + 13 211.836 5.000e-04 POINT (70.00000 70.00000) + 21 213.360 5.000e-04 POINT (30.00000 40.00000) >>> wn_gis.to_crs('EPSG:3857') >>> print(wn_gis.junctions.head()) - node_type elevation initial_quality geometry - name - 10 Junction 216.408 5.000e-04 POINT (2226389.816 11068715.659) - 11 Junction 216.408 5.000e-04 POINT (3339584.724 11068715.659) - 12 Junction 213.360 5.000e-04 POINT (5565974.540 11068715.659) - 13 Junction 211.836 5.000e-04 POINT (7792364.356 11068715.659) - 21 Junction 213.360 5.000e-04 POINT (3339584.724 4865942.280) + elevation initial_quality geometry + name + 10 216.408 5.000e-04 POINT (2226389.816 11068715.659) + 11 216.408 5.000e-04 POINT (3339584.724 11068715.659) + 12 213.360 5.000e-04 POINT (5565974.540 11068715.659) + 13 211.836 5.000e-04 POINT (7792364.356 11068715.659) + 21 213.360 5.000e-04 POINT (3339584.724 4865942.280) Snap point geometries to the nearest point or line ---------------------------------------------------- diff --git a/documentation/model_io.rst b/documentation/model_io.rst index f66f8fd53..bcd5b7636 100644 --- a/documentation/model_io.rst +++ b/documentation/model_io.rst @@ -206,27 +206,29 @@ GeoJSON files GeoJSON files are commonly used to store geographic data structures. More information on GeoJSON files can be found at https://geojson.org. -To use GeoJSON files in WNTR, a set of valid base column names are required. -Valid base GeoJSON column names can be obtained using the -:class:`~wntr.network.io.valid_gis_names` function. -The following example returns valid base GeoJSON column names for junctions. +When reading GeoJSON files into WNTR, only a set of valid column names can be used. +Valid GeoJSON column names can be obtained using the +:class:`~wntr.network.io.valid_gis_names` function. By default, the function +returns all column names, both required and optional. +The following example returns valid GeoJSON column names for junctions. .. doctest:: :skipif: gpd is None >>> geojson_column_names = wntr.network.io.valid_gis_names() >>> print(geojson_column_names['junctions']) - ['name', 'elevation', 'coordinates', 'emitter_coefficient', 'initial_quality', 'minimum_pressure', 'required_pressure', 'pressure_exponent', 'tag'] + ['name', 'elevation', 'geometry', 'emitter_coefficient', 'initial_quality', 'minimum_pressure', 'required_pressure', 'pressure_exponent', 'tag'] -A minimal list of valid column names can also be obtained by setting ``complete_list`` to False. -Column names that are optional (i.e., ``initial_quality``) and not included in the GeoJSON file are defined using default values. +A minimal list of required column names can also be obtained by setting ``complete_list`` to False. +Column names that are optional (i.e., ``initial_quality``) and not included in the GeoJSON file are +defined using default values. .. doctest:: :skipif: gpd is None >>> geojson_column_names = wntr.network.io.valid_gis_names(complete_list=False) >>> print(geojson_column_names['junctions']) - ['name', 'elevation', 'coordinates'] + ['name', 'elevation', 'geometry'] Note that GeoJSON files can contain additional custom column names that are assigned to WaterNetworkModel objects. @@ -253,7 +255,7 @@ Note that patterns, curves, sources, controls, and options are not stored in the The :class:`~wntr.network.io.read_geojson` function creates a WaterNetworkModel from a dictionary of GeoJSON files. -Valid base column names and additional custom attributes are added to the model. +Valid column names and additional custom attributes are added to the model. The function can also be used to append information from GeoJSON files into an existing WaterNetworkModel. .. doctest:: @@ -300,20 +302,21 @@ To use Esri Shapefiles in WNTR, several formatting requirements are enforced: assumed that the first 10 characters of each attribute are unique. * To create WaterNetworkModel from Shapefiles, a set of valid field names are required. - Valid base Shapefiles field names can be obtained using the - :class:`~wntr.network.io.valid_gis_names` function. - For Shapefiles, the `truncate` input parameter should be set to 10 (characters). - The following example returns valid base Shapefile field names for junctions. - Note that attributes like ``base_demand`` are truncated to ``base_deman``. + Valid Shapefiles field names can be obtained using the + :class:`~wntr.network.io.valid_gis_names` function. By default, the function + returns all column names, both required and optional. + For Shapefiles, the `truncate_names` input parameter should be set to 10 (characters). + The following example returns valid Shapefile field names for junctions. + Note that attributes like ``minimum_pressure`` are truncated to ``minimum_pr``. .. doctest:: :skipif: gpd is None >>> shapefile_field_names = wntr.network.io.valid_gis_names(truncate_names=10) >>> print(shapefile_field_names['junctions']) - ['name', 'elevation', 'coordinate', 'emitter_co', 'initial_qu', 'minimum_pr', 'required_p', 'pressure_e', 'tag'] + ['name', 'elevation', 'geometry', 'emitter_co', 'initial_qu', 'minimum_pr', 'required_p', 'pressure_e', 'tag'] - A minimal list of valid field names can also be obtained by setting ``complete_list`` to False. + A minimal list of required field names can also be obtained by setting ``complete_list`` to False. Field names that are optional (i.e., ``initial_quality``) and not included in the Shapefile are defined using default values. .. doctest:: @@ -322,7 +325,7 @@ To use Esri Shapefiles in WNTR, several formatting requirements are enforced: >>> shapefile_field_names = wntr.network.io.valid_gis_names(complete_list=False, ... truncate_names=10) >>> print(shapefile_field_names['junctions']) - ['name', 'elevation', 'coordinate'] + ['name', 'elevation', 'geometry'] * Shapefiles can contain additional custom field names that are assigned to WaterNetworkModel objects. @@ -349,7 +352,7 @@ Note that patterns, curves, sources, controls, and options are not stored in the The :class:`~wntr.network.io.read_shapefile` function creates a WaterNetworkModel from a dictionary of Shapefile directories. -Valid base field names and additional custom field names are added to the model. +Valid field names and additional custom field names are added to the model. The function can also be used to append information from Shapefiles into an existing WaterNetworkModel. .. doctest:: diff --git a/wntr/gis/network.py b/wntr/gis/network.py index d9c4e9d6b..dd1a633e7 100644 --- a/wntr/gis/network.py +++ b/wntr/gis/network.py @@ -99,14 +99,21 @@ def _create_gis(self, wn, crs: str = None, pumps_as_points: bool = False, Represent valves as points (True) or lines (False), by default False """ - def _extract_geodataframe(df, crs=None, links_as_points=False): - # Drop any column with all NaN + def _extract_geodataframe(df, crs=None, valid_base_names=None, + links_as_points=False): + if valid_base_names is None: + valid_base_names = [] + + # Drop any column with all NaN, this removes excess attributes + # Valid base attributes that have all None values are added back + # at the end of this routine df = df.loc[:, ~df.isna().all()] + # Define geom and drop node_type/link_type if df.shape[0] > 0: - # Define geom if 'node_type' in df.columns: geom = [Point((x,y)) for x,y in df['coordinates']] + del df['node_type'] elif 'link_type' in df.columns: geom = [] for link_name in df['name']: @@ -120,16 +127,25 @@ def _extract_geodataframe(df, crs=None, links_as_points=False): ls.append(v) ls.append(link.end_node.coordinates) geom.append(LineString(ls)) + del df['link_type'] - # Drop column if not a str, float, int, or bool + # Drop column if not a str, float, int, or bool (or np.bool_) + # This drops columns like coordinates, vertices # This could be extended to keep additional data type (list, # tuple, network elements like Patterns, Curves) drop_cols = [] for col in df.columns: - if not isinstance(df.iloc[0][col], (str, float, int, bool)): + # Added np.bool_ to the following check + # Returned by df.to_dict('records') for some network models + if not isinstance(df.iloc[0][col], (str, float, int, bool, np.bool_)): drop_cols.append(col) df = df.drop(columns=drop_cols) + # Add back in valid base attributes that had all None values + cols = list(set(valid_base_names) - set(df.columns)) + if len(cols) > 0: + df[cols] = None + # Set index if len(df) > 0: df.set_index('name', inplace=True) @@ -137,7 +153,7 @@ def _extract_geodataframe(df, crs=None, links_as_points=False): df = gpd.GeoDataFrame(df, crs=crs, geometry=geom) else: df = gpd.GeoDataFrame() - + return df # Convert the WaterNetworkModel to a dictionary @@ -146,29 +162,31 @@ def _extract_geodataframe(df, crs=None, links_as_points=False): df_nodes = pd.DataFrame(wn_dict['nodes']) df_links = pd.DataFrame(wn_dict['links']) + valid_base_names = self._valid_names(complete_list=False, truncate_names=None) + # Junctions df = df_nodes[df_nodes['node_type'] == 'Junction'] - self.junctions = _extract_geodataframe(df, crs) + self.junctions = _extract_geodataframe(df, crs, valid_base_names['junctions']) # Tanks df = df_nodes[df_nodes['node_type'] == 'Tank'] - self.tanks = _extract_geodataframe(df, crs) + self.tanks = _extract_geodataframe(df, crs, valid_base_names['tanks']) # Reservoirs df = df_nodes[df_nodes['node_type'] == 'Reservoir'] - self.reservoirs = _extract_geodataframe(df, crs) + self.reservoirs = _extract_geodataframe(df, crs, valid_base_names['reservoirs']) # Pipes df = df_links[df_links['link_type'] == 'Pipe'] - self.pipes = _extract_geodataframe(df, crs, False) + self.pipes = _extract_geodataframe(df, crs, valid_base_names['pipes'], False) # Pumps df = df_links[df_links['link_type'] == 'Pump'] - self.pumps = _extract_geodataframe(df, crs, pumps_as_points) + self.pumps = _extract_geodataframe(df, crs, valid_base_names['pumps'], pumps_as_points) # Valves df = df_links[df_links['link_type'] == 'Valve'] - self.valves = _extract_geodataframe(df, crs, valves_as_points) + self.valves = _extract_geodataframe(df, crs, valid_base_names['valves'], valves_as_points) def _create_wn(self, append=None): """ @@ -187,22 +205,32 @@ def _create_wn(self, append=None): wn_dict['nodes'] = [] wn_dict['links'] = [] - for element in [self.junctions, self.tanks, self.reservoirs]: + # Modifications to create a WaterNetworkModel from a dict + # Reset index + # Create coordinates/vertices from geometry + # Add node_type/link_type + for node_type, element in [('Junction', self.junctions), + ('Tank', self.tanks), + ('Reservoir', self.reservoirs)]: if element.shape[0] > 0: assert (element['geometry'].geom_type).isin(['Point']).all() df = element.reset_index(names="name") df.rename(columns={'geometry':'coordinates'}, inplace=True) df['coordinates'] = [[x,y] for x,y in zip(df['coordinates'].x, df['coordinates'].y)] + df['node_type'] = node_type wn_dict['nodes'].extend(df.to_dict('records')) - for element in [self.pipes, self.pumps, self.valves]: + for link_type, element in [('Pipe', self.pipes), + ('Pump', self.pumps), + ('Valve', self.valves)]: if element.shape[0] > 0: assert 'start_node_name' in element.columns assert 'end_node_name' in element.columns df = element.reset_index(names="name") df['vertices'] = df.apply(lambda row: list(row.geometry.coords)[1:-1], axis=1) df.drop(columns=['geometry'], inplace=True) + df['link_type'] = link_type wn_dict['links'].extend(df.to_dict('records')) # Create WaterNetworkModel from dictionary @@ -470,6 +498,17 @@ def _valid_names(self, complete_list=True, truncate_names=None): if truncate_names is not None and truncate_names > 0: for element, attributes in valid_names.items(): valid_names[element] = [attribute[:truncate_names] for attribute in attributes] + + for key, vals in valid_names.items(): + # Remove coordinates and vertices (not used to create GeoDataFrame geometry) + if 'coordinates' in valid_names[key]: + valid_names[key].remove('coordinates') + if 'vertices' in valid_names[key]: + valid_names[key].remove('vertices') + + # Add geometry + if 'geometry' not in valid_names[key]: + valid_names[key].append('geometry') return valid_names diff --git a/wntr/network/elements.py b/wntr/network/elements.py index 28e3e7fe8..24453ccd9 100644 --- a/wntr/network/elements.py +++ b/wntr/network/elements.py @@ -394,7 +394,7 @@ class Tank(Node): "min_level", "max_level", "diameter", - "min_vol" + "min_vol", "vol_curve_name", "overflow", "coordinates"] @@ -1041,7 +1041,7 @@ class Pump(Link): "end_node_name", "pump_type", "pump_curve_name", - "power" + "power", "base_speed", "speed_pattern_name", "initial_status"] diff --git a/wntr/network/io.py b/wntr/network/io.py index e4ceac8de..a299e459a 100644 --- a/wntr/network/io.py +++ b/wntr/network/io.py @@ -644,12 +644,13 @@ def valid_gis_names(complete_list=True, truncate_names=None): Valid column/field names for GeoJSON or Shapefiles Note that Shapefile field names are truncated to 10 characters - (set truncate=10) + (set truncate_names=10) Parameters ---------- complete_list : bool - Include a complete list of column/field names (beyond basic attributes) + When true, returns both optional and required column/field names. + When false, only returns required column/field names. truncate_names : None or int Truncate column/field names to specified number of characters, set truncate=10 for Shapefiles. None indicates no truncation. diff --git a/wntr/tests/test_gis.py b/wntr/tests/test_gis.py index 0872e179f..3018c80ae 100644 --- a/wntr/tests/test_gis.py +++ b/wntr/tests/test_gis.py @@ -132,12 +132,12 @@ def test_wn_to_gis(self): #assert self.gis_data.valves.shape[0] == self.wn.num_valves # Check minimal set of attributes - assert set(['node_type', 'elevation', 'geometry']).issubset(self.gis_data.junctions.columns) - assert set(['node_type', 'elevation', 'geometry']).issubset(self.gis_data.tanks.columns) - assert set(['node_type', 'geometry']).issubset(self.gis_data.reservoirs.columns) - assert set(['link_type', 'start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.pipes.columns) - assert set(['link_type', 'start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.pumps.columns) - #assert set(['link_type', 'start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.valves.columns) # Net1 has no valves + assert set(['elevation', 'geometry']).issubset(self.gis_data.junctions.columns) + assert set(['elevation', 'geometry']).issubset(self.gis_data.tanks.columns) + assert set(['geometry']).issubset(self.gis_data.reservoirs.columns) + assert set(['start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.pipes.columns) + assert set(['start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.pumps.columns) + #assert set(['start_node_name', 'end_node_name', 'geometry']).issubset(self.gis_data.valves.columns) # Net1 has no valves def test_gis_to_wn(self): @@ -265,15 +265,17 @@ def test_set_crs_to_crs(self): def test_add_attributes_and_write(self): - self.gis_data.add_node_attributes(self.results.node['pressure'].loc[3600,:], 'Pressure_1hr') - self.gis_data.add_link_attributes(self.results.link['flowrate'].loc[3600,:], 'Flowrate_1hr') + gis_data = self.wn.to_gis() + + gis_data.add_node_attributes(self.results.node['pressure'].loc[3600,:], 'Pressure_1hr') + gis_data.add_link_attributes(self.results.link['flowrate'].loc[3600,:], 'Flowrate_1hr') - assert 'Pressure_1hr' in self.gis_data.junctions.columns - assert 'Pressure_1hr' in self.gis_data.tanks.columns - assert 'Pressure_1hr' in self.gis_data.reservoirs.columns - assert 'Flowrate_1hr' in self.gis_data.pipes.columns - assert 'Flowrate_1hr' in self.gis_data.pumps.columns - assert 'Flowrate_1hr' not in self.gis_data.valves.columns # Net1 has no valves + assert 'Pressure_1hr' in gis_data.junctions.columns + assert 'Pressure_1hr' in gis_data.tanks.columns + assert 'Pressure_1hr' in gis_data.reservoirs.columns + assert 'Flowrate_1hr' in gis_data.pipes.columns + assert 'Flowrate_1hr' in gis_data.pumps.columns + assert 'Flowrate_1hr' not in gis_data.valves.columns # Net1 has no valves def test_write_geojson(self): prefix = 'temp_Net1' diff --git a/wntr/tests/test_network.py b/wntr/tests/test_network.py index 5e30e1204..f6428ecf2 100644 --- a/wntr/tests/test_network.py +++ b/wntr/tests/test_network.py @@ -1087,6 +1087,30 @@ def test_shapefile_roundtrip(self): files['valves'] = 'temp_valves' B = self.wntr.network.read_shapefile(files) assert(wn._compare(B, level=0)) - + + def test_valid_gis_names(self): + + required_names = wntr.network.io.valid_gis_names(complete_list=False, truncate_names=None) + valid_names = wntr.network.io.valid_gis_names(complete_list=True, truncate_names=None) + + wn = self.wntr.network.WaterNetworkModel(join(ex_datadir, "Net6.inp")) + gis_data = wn.to_gis() + + for component in required_names.keys(): + required_columns = required_names[component] + valid_columns = valid_names[component] + + data = getattr(gis_data, component) + data_columns = list(data.columns) + data_columns.append(data.index.name) + + # Check that all data columns are valid + assert len(set(data_columns)-set(valid_columns)) == 0 + # Check that all required columns are in the data + assert len(set(required_columns)-set(data_columns)) == 0 + # Assert that node_type and link_type are not in data columns + assert 'node_type' not in data_columns + assert 'link_type' not in data_columns + if __name__ == "__main__": unittest.main() From 03cf055fc1f33aff05182769375d5a179ab6c2ee Mon Sep 17 00:00:00 2001 From: Katherine Klise Date: Tue, 12 Nov 2024 09:23:03 -0800 Subject: [PATCH 5/5] Bug fix in roughness unit conversion when using D-W (#450) * Added D-W conversion for roughness coefficient for inp file i/o * bug fix in roughness conversion * Added a warning if the user changes headloss to/from D-W * Added additional tests for changing headloss formula * Updated test * commented out graphics in test --- wntr/epanet/io.py | 11 ++++- wntr/epanet/util.py | 4 +- wntr/network/options.py | 14 +++++++ wntr/tests/test_epanet_io.py | 79 ++++++++++++++++++++++++++++++++++-- 4 files changed, 101 insertions(+), 7 deletions(-) diff --git a/wntr/epanet/io.py b/wntr/epanet/io.py index fea4e3651..af56b8bed 100644 --- a/wntr/epanet/io.py +++ b/wntr/epanet/io.py @@ -674,6 +674,8 @@ def _write_tanks(self, f, wn, version=2.2): f.write('\n'.encode(sys_default_enc)) def _read_pipes(self): + darcy_weisbach = self.wn.options.hydraulic.headloss == "D-W" + for lnum, line in self.sections['[PIPES]']: line = line.split(';')[0] current = line.split() @@ -701,7 +703,8 @@ def _read_pipes(self): current[2], to_si(self.flow_units, float(current[3]), HydParam.Length), to_si(self.flow_units, float(current[4]), HydParam.PipeDiameter), - float(current[5]), + to_si(self.flow_units, float(current[5]), HydParam.RoughnessCoeff, + darcy_weisbach=darcy_weisbach), minor_loss, link_status, check_valve) @@ -711,6 +714,8 @@ def _read_pipes(self): raise ENValueError(211, str(e.args[0]), line_num=lnum) from e def _write_pipes(self, f, wn): + darcy_weisbach = wn.options.hydraulic.headloss == "D-W" + f.write('[PIPES]\n'.encode(sys_default_enc)) f.write(_PIPE_LABEL.format(';ID', 'Node1', 'Node2', 'Length', 'Diameter', 'Roughness', 'Minor Loss', 'Status').encode(sys_default_enc)) @@ -723,7 +728,9 @@ def _write_pipes(self, f, wn): 'node2': pipe.end_node_name, 'len': from_si(self.flow_units, pipe.length, HydParam.Length), 'diam': from_si(self.flow_units, pipe.diameter, HydParam.PipeDiameter), - 'rough': pipe.roughness, + 'rough': from_si(self.flow_units, pipe.roughness, + HydParam.RoughnessCoeff, + darcy_weisbach=darcy_weisbach), 'mloss': pipe.minor_loss, 'status': str(pipe.initial_status), 'com': ';'} diff --git a/wntr/epanet/util.py b/wntr/epanet/util.py index c06528573..c09e51832 100644 --- a/wntr/epanet/util.py +++ b/wntr/epanet/util.py @@ -571,7 +571,7 @@ def _to_si(self, flow_units, data, darcy_weisbach=False): elif self in [HydParam.RoughnessCoeff] and darcy_weisbach: if flow_units.is_traditional: - data = data * (1000.0 * 0.3048) # 1e-3 ft to m + data = data * (0.001 * 0.3048) # 1e-3 ft to m elif flow_units.is_metric: data = data * 0.001 # mm to m @@ -668,7 +668,7 @@ def _from_si(self, flow_units, data, darcy_weisbach=False): elif self in [HydParam.RoughnessCoeff] and darcy_weisbach: if flow_units.is_traditional: - data = data / (1000.0 * 0.3048) # 1e-3 ft from m + data = data / (0.001 * 0.3048) # 1e-3 ft from m elif flow_units.is_metric: data = data / 0.001 # mm from m diff --git a/wntr/network/options.py b/wntr/network/options.py index 3d1d086a5..120cb39dc 100644 --- a/wntr/network/options.py +++ b/wntr/network/options.py @@ -10,6 +10,7 @@ """ import re import logging +import warnings import copy logger = logging.getLogger(__name__) @@ -406,6 +407,19 @@ def __setattr__(self, name, value): value = str.upper(value) if value not in ['H-W', 'D-W', 'C-M']: raise ValueError('headloss must be one of "H-W", "D-W", or "C-M"') + # If headloss is changed from ['H-W', 'C-M'] to/from 'D-W', print + # a warning, the units of the roughness coefficient cannot be + # converted from unitless to length (0.001 ft or mm) + try: + orig_value = self.__dict__[name] + except: + orig_value = None + if orig_value is not None: + if (orig_value in ['H-W', 'C-M'] and value == 'D-W') or \ + (value in ['H-W', 'C-M'] and orig_value == 'D-W'): + warnings.warn('Changing the headloss formula from ' + + orig_value + ' to ' + value + + ' will not change the units of the roughness coefficient.') elif name == 'hydraulics': if value is not None: value = str.upper(value) diff --git a/wntr/tests/test_epanet_io.py b/wntr/tests/test_epanet_io.py index 38d611734..8414ca127 100644 --- a/wntr/tests/test_epanet_io.py +++ b/wntr/tests/test_epanet_io.py @@ -625,7 +625,9 @@ def setUpClass(self): def tearDownClass(self): pass - def test_link_flowrate_units_convert(self): + def test_units_convert(self): + # Compares Net3 EpanetSimulator flowrate results using INP files saved + # using GPM and CMH units for link_name, link in self.wn.links(): for t in self.results2.link["flowrate"].index: self.assertLessEqual( @@ -636,7 +638,7 @@ def test_link_flowrate_units_convert(self): 0.00001, ) - def test_link_headloss_units_convert(self): + def test_link_headloss_units(self): # headloss = per unit length for pipes and CVs pipe_name = '123' @@ -665,7 +667,78 @@ def test_link_headloss_units_convert(self): ), 0.00001, ) - + def test_headloss_formula_roughness_units(self): + """ + See Table 3.2 Roughness Coefficients, this test uses values for Plastic Pipes + https://epanet22.readthedocs.io/en/latest/3_network_model.html?highlight=roughness#id3 + """ + pressure = {} + for headloss in ['H-W', 'C-M', 'D-W']: + for units in ['GPM', 'LPS']: + file_prefix = 'temp_'+headloss+'_'+units + + inp_file = join(ex_datadir, "Net3.inp") + wn = self.wntr.network.WaterNetworkModel(inp_file) + wn.options.hydraulic.demand_model = 'DDA' + wn.options.hydraulic.inpfile_units = units + wn.options.time.duration = 0 # steady state + + wn.options.hydraulic.headloss = headloss + if headloss == 'H-W': + roughness = 145 # unitless + elif headloss == 'D-W': + roughness = 0.005 # 0.001*ft + roughness = (roughness*0.001)*0.3048 # 1.524e-6 m == 0.001524 mm + else: # C-M + roughness = 0.011 # unitless + + for name, pipe in wn.pipes(): + pipe.roughness = roughness + + sim = self.wntr.sim.EpanetSimulator(wn) + results = sim.run_sim(file_prefix=file_prefix) + temp = results.node['pressure'].loc[0,wn.junction_name_list] + pressure[headloss+', '+units] = temp + + #import matplotlib.pylab as plt + #import pandas as pd + #plt.figure() + #pd.DataFrame(pressure).plot.bar() + #plt.legend() #loc='lower right') + + ## Compares all results to H-W, GPM + threshold = 0.55 + for key in pressure.keys(): + MAE = (pressure['H-W, GPM'] - pressure[key]).abs().mean() + print(key, MAE) + self.assertLessEqual(MAE, threshold) # m + + ## Changing the headloss formula from H-W to D-W, without changing roughness, + # will result in errors + inp_file = join(ex_datadir, "Net3.inp") + wn = self.wntr.network.WaterNetworkModel(inp_file) + assert wn.options.hydraulic.headloss == 'H-W' + assert wn.options.hydraulic.inpfile_units == 'GPM' + pressure_hw = pressure['H-W, GPM'] + # Change to D-W without changing roughness + wn.options.hydraulic.headloss = 'D-W' + wn.options.hydraulic.demand_model = 'DDA' + wn.options.hydraulic.inpfile_units = units + wn.options.time.duration = 0 # steady state + sim = self.wntr.sim.EpanetSimulator(wn) + results = sim.run_sim() + pressure_dw = results.node['pressure'].loc[0,wn.junction_name_list] + # Compare results + MAE = (pressure_hw - pressure_dw).abs().mean() + with self.assertRaises(AssertionError): + self.assertLessEqual(MAE, threshold) # m + + ## D-W is not supported by the WNTRSimulator + sim = self.wntr.sim.WNTRSimulator(wn) + with self.assertRaises(NotImplementedError): + results = sim.run_sim() + + if __name__ == "__main__": unittest.main()