maxflow/get_layerBudget.py at code_refactoring · KWB-R/maxflow · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
# -*- coding: utf-8 -*-
"""
Created on Mon Jan 16 10:58:55 2017

@author: mrustl
"""
import os
import numpy as np
import pandas as pd
import flopy.utils.binaryfile as bf


def get_layerbudget(modelname,
                    nper,
                    perlen,
                    nlay,
                    model_ws = '.',
                    debug = True
                    ):

    perlen = np.array(perlen, ndmin=1, copy=False)

    bud_agg = pd.DataFrame(columns=['stress_period',
                              'time_step',
                              'layer',
                              'STORAGE_IN',
                              'STORAGE_OUT',
                              'CONSTANT_HEAD_IN',
                              'CONSTANT_HEAD_OUT',
                              'FLOW_RIGHT_FACE',
                              'FLOW_FRONT_FACE',
                              'FLOW_LOWER_FACE'])
    cbb = bf.CellBudgetFile(os.path.join(model_ws, modelname+'.cbc'))
    for stress_period in np.arange(0,nper).astype('int'):
        for time_step in np.arange(0,perlen[stress_period]).astype('int'):
            bud = cbb.get_data(kstpkper = (time_step,stress_period),
                       full3D=True)
            for layer in np.arange(0,nlay).astype('int'):
                if debug: print("Stress period: " + str(stress_period + 1) + ", Time step: " + str(time_step + 1) + ", Layer: " + str(layer + 1))
                tmp = pd.DataFrame([[stress_period,
                         time_step + 1,
                         layer,
                         bud[0][layer][bud[0][layer]>0].sum(),
                         bud[0][layer][bud[0][layer]<0].sum(),
                         bud[1][layer][bud[1][layer]>0].sum(),
                         bud[1][layer][bud[1][layer]<0].sum(),
                         bud[2][layer].sum(),
                         bud[3][layer].sum(),
                         bud[4][layer].sum()
                         ]],
                     columns=['stress_period',
                              'time_step',
                              'layer',
                              'STORAGE_IN',
                              'STORAGE_OUT',
                              'CONSTANT_HEAD_IN',
                              'CONSTANT_HEAD_OUT',
                              'FLOW_RIGHT_FACE',
                              'FLOW_FRONT_FACE',
                              'FLOW_LOWER_FACE'])
                bud_agg = bud_agg.append(tmp,ignore_index=True)
    bud_agg.loc[:,['CONSTANT_HEAD_IN']] = bud_agg['CONSTANT_HEAD_IN'].as_matrix().astype("float32")
    bud_agg.loc[:,['CONSTANT_HEAD_OUT']] = bud_agg['CONSTANT_HEAD_OUT'].as_matrix().astype("float32")
    bud_agg.loc[np.isnan(bud_agg['CONSTANT_HEAD_IN']),['CONSTANT_HEAD_IN']] = 0
    bud_agg.loc[np.isnan(bud_agg['CONSTANT_HEAD_OUT']),['CONSTANT_HEAD_OUT']] = 0
    return(bud_agg)