From eef5d333beda3a69e6cdd8d0cfd257f46671e6ed Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Tue, 26 Nov 2024 14:23:01 -0500 Subject: [PATCH] feat(binaryfile): add head/budget file reversal script --- .../{binaryfile.py => binaryfile/__init__.py} | 96 +++++++++--------- flopy/utils/binaryfile/reverse.py | 31 ++++++ flopy/utils/util_array.py | 2 +- freyberg-rev.hds | Bin 0 -> 6452 bytes 4 files changed, 80 insertions(+), 49 deletions(-) rename flopy/utils/{binaryfile.py => binaryfile/__init__.py} (97%) create mode 100644 flopy/utils/binaryfile/reverse.py create mode 100644 freyberg-rev.hds diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile/__init__.py similarity index 97% rename from flopy/utils/binaryfile.py rename to flopy/utils/binaryfile/__init__.py index 4aa31260d..82edb8830 100644 --- a/flopy/utils/binaryfile.py +++ b/flopy/utils/binaryfile/__init__.py @@ -19,8 +19,8 @@ import numpy as np import pandas as pd -from ..utils.datafile import Header, LayerFile -from .gridutil import get_lni +from ..datafile import Header, LayerFile +from ..gridutil import get_lni HEAD_TEXT = " HEAD" @@ -663,38 +663,28 @@ def get_max_kper_kstp_tsim(): kstp[header["kper"]] = 0 return kper, kstp, tsim - # get max period and time from the head file maxkper, maxkstp, maxtsim = get_max_kper_kstp_tsim() - # if we have tdis, get max period number and simulation time from it - tdis_maxkper, tdis_maxtsim = None, None - if self.tdis is not None: - pd = self.tdis.perioddata.get_data() - if any(pd): - tdis_maxkper = len(pd) - 1 - tdis_maxtsim = sum([p[0] for p in pd]) - # if we have both, check them against each other - if tdis_maxkper is not None: - assert maxkper == tdis_maxkper, ( - f"Max stress period in binary head file ({maxkper}) != " - f"max stress period in provided tdis ({tdis_maxkper})" - ) - assert maxtsim == tdis_maxtsim, ( - f"Max simulation time in binary head file ({maxtsim}) != " - f"max simulation time in provided tdis ({tdis_maxtsim})" - ) + prev_kper = None + perlen = None def reverse_header(header): """Reverse period, step and time fields in the record header""" + nonlocal prev_kper + nonlocal perlen + # reverse kstp and kper headers kstp = header["kstp"] - 1 kper = header["kper"] - 1 header["kstp"] = maxkstp[kper] - kstp + 1 header["kper"] = maxkper - kper + 1 + if kper != prev_kper: + perlen = header["pertim"] + prev_kper = kper + # reverse totim and pertim headers header["totim"] = maxtsim - header["totim"] - perlen = pd[kper][0] header["pertim"] = perlen - header["pertim"] return header @@ -1022,7 +1012,6 @@ def __init__( self.paknamlist_from = [] self.paknamlist_to = [] self.compact = True # compact budget file flag - self.dis = None self.modelgrid = None if "model" in kwargs.keys(): @@ -2237,24 +2226,46 @@ def reverse(self, filename: Optional[os.PathLike] = None): ] ) - # make sure we have tdis - if self.tdis is None or not any(self.tdis.perioddata.get_data()): - raise ValueError("tdis must be known to reverse a cell budget file") + nrecords = len(self) + target = filename + + def get_max_kper_kstp_tsim(): + header = self.recordarray[-1] + kper = header["kper"] - 1 + tsim = header["totim"] + kstp = {0: 0} + for i in range(len(self) - 1, -1, -1): + header = self.recordarray[i] + if header["kper"] in kstp and header["kstp"] > kstp[header["kper"]]: + kstp[header["kper"]] += 1 + else: + kstp[header["kper"]] = 0 + return kper, kstp, tsim + + maxkper, maxkstp, maxtsim = get_max_kper_kstp_tsim() + prev_kper = None + perlen = None - # extract perioddata - pd = self.tdis.perioddata.get_data() + def reverse_header(header): + """Reverse period, step and time fields in the record header""" - # get maximum period number and total simulation time - nper = len(pd) - kpermx = nper - 1 - tsimtotal = 0.0 - for tpd in pd: - tsimtotal += tpd[0] + nonlocal prev_kper + nonlocal perlen - # get number of records - nrecords = len(self) + # reverse kstp and kper headers + kstp = header["kstp"] - 1 + kper = header["kper"] - 1 + header["kstp"] = maxkstp[kper] - kstp + 1 + header["kper"] = maxkper - kper + 1 - target = filename + if kper != prev_kper: + perlen = header["pertim"] + prev_kper = kper + + # reverse totim and pertim headers + header["totim"] = maxtsim - header["totim"] + header["pertim"] = perlen - header["pertim"] + return header # if rewriting the same file, write # temp file then copy it into place @@ -2269,18 +2280,7 @@ def reverse(self, filename: Optional[os.PathLike] = None): for idx in range(nrecords - 1, -1, -1): # load header array header = self.recordarray[idx] - - # reverse kstp and kper in the header array - (kstp, kper) = (header["kstp"] - 1, header["kper"] - 1) - kstpmx = pd[kper][1] - 1 - kstpb = kstpmx - kstp - kperb = kpermx - kper - (header["kstp"], header["kper"]) = (kstpb + 1, kperb + 1) - - # reverse totim and pertim in the header array - header["totim"] = tsimtotal - header["totim"] - perlen = pd[kper][0] - header["pertim"] = perlen - header["pertim"] + header = reverse_header(header) # Write main header information to backward budget file h = header[ diff --git a/flopy/utils/binaryfile/reverse.py b/flopy/utils/binaryfile/reverse.py new file mode 100644 index 000000000..f69e00f5d --- /dev/null +++ b/flopy/utils/binaryfile/reverse.py @@ -0,0 +1,31 @@ +import argparse +from pathlib import Path + +from flopy.utils.binaryfile import CellBudgetFile, HeadFile + +if __name__ == "__main__": + """Reverse head or budget files.""" + + parser = argparse.ArgumentParser(description="Reverse head or budget files.") + parser.add_argument( + "--infile", + "-i", + type=str, + help="Input file.", + ) + parser.add_argument( + "--outfile", + "-o", + type=str, + help="Output file.", + ) + args = parser.parse_args() + infile = Path(args.infile) + outfile = Path(args.outfile) + suffix = infile.suffix.lower() + if suffix in [".hds", ".hed"]: + HeadFile(infile).reverse(outfile) + elif suffix in [".bud", ".cbc"]: + CellBudgetFile(infile).reverse(outfile) + else: + raise ValueError(f"Unrecognized file suffix: {suffix}") diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py index 41c0e4410..6d117b840 100644 --- a/flopy/utils/util_array.py +++ b/flopy/utils/util_array.py @@ -14,8 +14,8 @@ import numpy as np from ..datbase import DataInterface, DataType -from ..utils.binaryfile import BinaryHeader from ..utils.flopy_io import line_parse +from .binaryfile import BinaryHeader class ArrayFormat: diff --git a/freyberg-rev.hds b/freyberg-rev.hds new file mode 100644 index 0000000000000000000000000000000000000000..71c1473a482466d2315e46fe7def2bf784db721d GIT binary patch literal 6452 zcmbW63sg;O*T=i4(}f};Qc|ehuD#dVYwcTm$Dc$=hbSsiDcx_%(S<^kh$LMdktlZ& z9STWAE=ffcsa&EGB9znRTkpj?zVW`}4Da_nW2`a9*kg^ce|!GtoX?zd%gV^e$TH8r zAHDx|tSo2F9>cs~#5{j{SD>Vo-#1HwSBmZN!G;pFD0#U`JItV$bLVx;ZZmA!?B?%b zXog1(?$i8^h>>Qb#o9YbjLUw!t7Fqdh%jj?s~r^LeNknoj;9b-c@_^7E(?$nux8`d zQ37b{-7c+O#>d&_ZMS_-(y*LXq5mH8=fxJ6GK#LDKnTN!lL(9CkDR>agcnMdU=r$UJYtRBaaI3sV%v!YJiVwen;8? zOFo>xWX-YPN#j9~cv#wN3LCi{o^5>;yp~4vhK!)$({gpX(gYd}NuC;7`ZVSn2Kr_j z(x}cV+v8d&tFMq8u|kd>3v!+uRIV{ySnn6p68PX3-ln zeKEcjEWcHhBtnpyXQWf75ZPNDo(zhFaLAkU@z`bo?CTdr4b|~+#gY?yX)GUhzj$uY z@}yC#dZ&KvDGECCQ!OhVP>4y;P|5G4(6|4+-l={HVe_^O+9=ToKewiEts0FjUYojt z)=Mz{(6;5p&Jx_YBC<3dEkV3S<<5geX0XrC%l2MihTERo7ggO5L&f8lUzbi1qxM^8 z#ltNks8oMlH|eerT~9u@y2uN$PQCt+uZ;jcSNon1?Bipcbo9N}b{cA5KiT=xjIO)o zOg+LV1Zx}I@7qgZ4(H2;gi91GZ+u_zwT{BH_%W*6S}E))bJz83r%>lQqk?}#g5;#7 z``&D1I91}OM2RE_dR$@jvB3p2=%lC-vy-zkKjsPM2CQ;0m3 zyJgP`3MMnCG)@c!{a@SGRqvy4$lGX!?pX?NN39w0C&N8jNBu@*y#!9xbEnKNl0bKL z#0r+b1bXwsI#Y*Bz>5qo(n~ahQR2$X&7ES{Tn?B-u*LA5;i;PSs|a`N-jwe35rSn| zs68T90MGC0{wcTlm>llCFy4R* zTY{qa(CkgVZ&pKN?)v_3R!lz5S+8t&X)6U~9hJ;aR8K-rkgolfYJsZBh;URct z-MNJE6wXFnPrGVHA+D!2Q*%BAXSqa+R&yz``4_DHrbw~K`F4l>2ML~RUUlYWssuI5 z*Sm@ynt`y3%`%oVgJ_0biHD9DqvyVMsy!ow${6Yb7{5^y!r!PA@^Q=GXk^1EK6chu zXKM!2=rn7uOngG2nO0vE!T8BpRilR=$HN#_uTirW507+i>rD&eAwj&xWBM%~E<|lg z#&;g-pH)0~Fp9#uyxn&!W26|KLix1Bic1~4%KFX=w7oj;cCN3!{&K^c%RaQ8EN;s3B$E$TNf8k-B zd_(b{b9r!(m6iWx9RWwn+C?heBwWVZ9DNbPgMy25U_RsD?EDmNWfKp6XTtfTN~AEn znZ8o_s~5VrRuDrx%v9GX<{lssI`Q_)3FCMemUW30+eE@(;P{9A{UrP&Bd6?V zF>!dV>FyyH9+tC%u5uXttv}1_yvgFB;Kmr&=?dluI&*nYC})nq6negHzZ9ZViqHJM z{haTjeVOIE0SY|w^EmQ81$oZ{lO_dGu-$$5f?hKZZJ#q9rO`a>`CO?ui}9Cx4BP5m zHi@R{*ag#@NhB6Dh;}k@c(H2T%>+vx5?5ZRH}vI!%5MqEi{ydVeg9D8D03{d$gG+_ z(j1g~-g(Or=FpR0qO?%i{4eA0hrWqUxl67WQK+>o*%p(?j_& zX`be$VhbKTw`g~<>PdX+(8SuHTU+ zDgK!8Sh=rJiU||+q8BpuhmL<`4~7{`jhscn!I^U*nyIhdn#)4hR`Nja)(=&8=fBiYD<|CV!y%B8h2vA%<>^Bu*v2 ztkiC2_@gxUo-35XmK9~5Q6NQz*WD3?Or1Fcvz+wt)sx?#pQ?%kkk+9He3>!mB@$q5c#PM7yPUYL| zd0-?3J$7nUC3D~Xm99TNZ=S?HsJo29*iS*boOLKH-nY0ilkwk%cY|qn2T64AO_%Xz zbTmyE3|`?#AmsPvc^WDNT&kP8&&3lMH2?O8{yPGV6@zOMX%e4L`Z+IKK|<^Mw%gp* zj2}zcP0yZ7Fl~*S#jHjNTqicMY92^nt4_X7tC0MZZ_99{ij8|HXlLZVKW{_9FEe+G z>`flVD(JqDVd`CUlG5vd8b-(a=1o4EIN-8Xlikv}@DF%jsp>#r^>u~5+$#jEhFPa7 zE0b`YZaGm~NWyJ;b@xVV5>dN#XW7R}@ItVJiBcM^iRS! z{Ld7~;r@penEOoWKwx7=4-fj~LigV%@(>yO;AGDbfvAF{hl^iv5uR;iVZh?za@O6{ z0eu4YhyN_ENG5Pjrm?i^8S^{N-E8HIAfa+{ZF8?aiB*p$CX(6`tWvVwAEhjTTeTi@PLGkh)h^4krr86?3aY7xxb z_Jgl$cD2=YWdRh`1Ujla`FL=KI~3AD<5Es(Srwbcz1*}1vP_*Gy<>QI*7_B4_uww1$}`j z?u{6?ea@>|-WKDMK+C}}MT|cRp584q7Gv<}h9TB^5ybcMhQdw@@m;xd?SvKq`Z6bJ zQR4+D9o9waiup*ib9|e=h(-=4O)e#r!jkH~n;)2dWOhb?#mI3a?#DzOvD(iC(!Ums zax=y5&gn5L$~YLcBgt*nOD@LQev8oi%m6}JMM&=ATK@3F=R3atF~;J z*5k@W?a_W#(nzKb?i12GtqFuKOxp2yC4mg_NpJs;BKWq%K8U?9!n#dsKCap;!nV;- z?Qbna7+3n?zSk2WqQ=U7iT4n~e(YO|rY8cty@|N@!vzSaFn4=s%E#ctZJ)(+X?#;_ z9o9U8LT~py(?O=6x^U1l_$E`|_SMy&Uig&*UQfwX2RRd1YixWcx^D`Hf<2b9MvPxs zZB_j}TnwZi^0DI)Xb?WjdODxL_HW%=gN}<3=#!qCmMDU9OXK*fP9ij$zw%8{6=7+@ zYqdRRg_yeY$ToMe5c$_?hwpeP!00%Mwq_3>D9qm-P)|d5vX;uI;WR2-{71AtgQnNpZa5lRdG=| zu-z|Sm%!4g)q`~c0?r2pN|psO@%6;bta_0MvKFIj!u3Q5uR5)5_FM=ND>XbIQ;0Xm z?!1;R5JD}qf_saZ!>f)|U-!DphmE?~^|%wvoc5bwXm&n@{i_pS%C2VmUM?2~0ury> zlAMi2j7~4FpLH6^)RjexA+ODbtjZp`Yl10cO0M*BFLAIzV}5aa1{d|EvF7XK3H1AhTnwD<^HQy2v;&)#~_7<(O?ZB)xm5WOqt`94<;{MOAj^S}_;7i&o6eX6pB{ zWlf=`1|poY&w3_)FT{mi{wK#AWBkSTi#xqkh~_z)jWS0I(Xm|L@^rZXGCdE=E=?35 z)VUzWe;gkbam&6GF#D$P&6K(ZGjFh!C)va>y8dq69z2Pu7eTwXcBfwBAfHZ9zSm`f zgkjNdXPOuzRH#JX6`Ej9=Eb4(3R7?n`)E2{=HQ8wom{^;7j((|=K7V~f87fk?|xNp zJ}1PKLr;lC?u@SQ-|`0(g_v`EqO3Qbb z?remf+mW-R-vs`-2kaKAbFhg#Vq2WgL8@EevOCW>|Kj^`e^oqH-O;2XgtM7PJ^z>h z9Ypp@tyuy@y2ML4%wFYs!lUPh1bqDNd3l0j3Nz=%?yp|Hp4p?IF5=o`9&T)?y?05D zK!5S}QNIhPB1Vw3RDXss5*i%V5n;wSwBTdh^jR!SR4gx3R$xP=Ljd9vw zS=ja}Be*=j1H)I4x#mDm`M4RoSi(%9hIf6KYF#ROVQPd+Y^ zH~kM?fAIXfb$@+C??XHd_jQ-uhrMI=U@|p9-u@I`^NY&nz2hN>4h%{4XXf0EiPO#V zNF3V84_g^apz`qai$hO2Aa6B)SMf6e`|6-ura}wW_6+v_nd=8%VD*L~&q`))xiqE!STcz@5$DR+%_Ol|yHZ#&NZ@AM zktxnL1b%sSF+^rA7x`Xw38RjgLWLWkuU*GROX_&el>!#5W36QZi`MJ$xe zd6RB8-T42_^@GprZ>NpG->nx@KeCv^^uyHO5BHyz^Dq04|2RzW4{Q8S%lfA^`mcHZ E2Q0sf0{{R3 literal 0 HcmV?d00001