-
Notifications
You must be signed in to change notification settings - Fork 0
/
adjust_topg
executable file
·99 lines (83 loc) · 3.57 KB
/
adjust_topg
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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#!/usr/bin/env python
import netCDF4 as nc
# import matplotlib.pylab as mpl
import sys
import flo_utils as fu
from argparse import ArgumentParser
import shutil
import numpy as np
def adjust_topo(filename, reference, state, opts):
dataset = nc.Dataset(filename, "r+")
reference = nc.Dataset(reference, "r")
state = nc.Dataset(state, "r+")
usurf_target = np.squeeze(reference.variables["usurf"][:])
reference_thickness = np.squeeze(reference.variables["thk"][:])
current_topg = np.squeeze(dataset.variables["topg"][:])
old_factors = np.squeeze(state.variables["factor"][:])
if opts["target"] == "smb":
current_smb_orig = np.squeeze(dataset.variables["climatic_mass_balance_original"][:])
current_smb = np.squeeze(dataset.variables["climatic_mass_balance"][:])
current_smb_correction = current_smb - current_smb_orig
# Wenn die Korrektur negativ ist, habe ich zu viel Eis, d.h. das Eis muss dicker werden, damit es schneller abfliesst, d.h. ich baue mir einen positiven Faktor.
# Deshalb - in dieser Gleichung.
factors = - current_smb_correction/100./910. # ice density...
else:
current_csurf = np.squeeze(dataset.variables["csurf"][:])
reference_csurf = np.squeeze(reference.variables["csurf"][:])
factors = - ((current_csurf - reference_csurf) / (current_csurf + reference_csurf)/1000.).filled(0.) # damping factor of 1000. arbitrary choice
# if fu.debug:
# mpl.figure()
# mpl.imshow(factors)
# mpl.title("error term")
# mpl.colorbar()
factors = factors + (factors==0)*1e-30
factors = (factors > 0 )* (1+factors) + (factors<0)/(1-factors)
# if fu.debug:
# mpl.figure()
# mpl.imshow(factors)
# mpl.title("raw factor change")
# mpl.colorbar()
factors = factors * old_factors
factors = (factors < 2) * (factors > .5) * factors + (factors >=2) * 2 + (factors <=0.5)*0.5
# if fu.debug:
# mpl.figure()
# mpl.imshow(factors)
# mpl.title("final factors")
# mpl.colorbar()
new_thk = reference_thickness * factors
# if fu.debug:
# mpl.figure()
# mpl.imshow(- (new_thk - reference_thickness))
# mpl.title("topg change")
# mpl.colorbar()
# mpl.show()
new_topg = usurf_target - new_thk
dataset.variables["topg"][0,:] = (new_topg*(new_topg <=usurf_target))+ (usurf_target*(new_topg >usurf_target))
dataset.variables["thk"][0,:] = (dataset.variables["thk"][0,:] - (dataset.variables["topg"][:] - current_topg))
dataset.variables["thk"][0,:] = (dataset.variables["thk"][:])*(dataset.variables["thk"][:]>0)
state.variables["factor"][0,:] = factors
def parse_args():
parser = ArgumentParser()
parser.description = "compare slopes of two variables from two files"
parser.add_argument("FILES", nargs=1)
parser.add_argument("-v", "--verbose",
help='''Be verbose''', action="store_true")
parser.add_argument("-r", "--reference",
help='''file with reference values''', default="reference.nc")
parser.add_argument("-s", "--state",
help='''file with reference values''', required = True)
parser.add_argument("-t", "--target",
help='''reference variable''', default = "csurf")
# parser.add_argument("-b", "--var_b",
# help='''variable b''', default="data")
options = parser.parse_args()
return options
def main():
options = parse_args()
if options.verbose:
print (dir(options))
print options.FILES
fu.debug=True
adjust_topo(options.FILES[0], options.reference, options.state, vars(options))
if __name__ == "__main__":
main()