Skip to content
Snippets Groups Projects
Commit ab62011a authored by Martin Schrön's avatar Martin Schrön
Browse files

New full example to vertically weight sm time serieses

parent 718a2e33
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/weighting/CRNS-weighting.py 0 → 100644
+ 165
− 0
View file @ ab62011a
# %%
import numpy as np
import pandas
import matplotlib.pyplot as plt
from Schroen2017hess import Wr, D86, Weight_d, get_footprint
# %%
# input values
class Profile:
__slots__ = [
# Input
"sm", # soil moisture values
"d", # depth values in cm
"bd", # bulk density
"r", # distance from the CRNS in m
# Calculated
"D86", # penetration depth
"w_r", # radial weight of this profile
"sm_wavg", # vertically weighted average sm
"data" # DataFrame
]
def __init__(self, sm, d, bd, r=1):
# Vector data
self.sm = np.array(sm)
self.d = np.array(d)
self.bd = np.array(bd)
self.data = None
self.update_data()
# Scalar values
self.r = r
self.D86 = np.nan
self.sm_wavg = np.nan
self.w_r = np.nan
def update_data(self):
if not self.data:
self.data = pandas.DataFrame()
self.data["sm"] = self.sm
self.data["d"] = self.d
self.data["bd"] = self.bd
if not "weight" in self.data.columns:
self.data["weight"] = np.nan
# %%
# Site
air_humidity = 5
air_pressure = 1013
profiles = [
Profile(
sm = [0.21, 0.42, 0.23, 0.19],
d = [ 5, 10, 20, 40],
bd = [ 1.5, 1.5, 1.5, 1.5],
r = 5
),
Profile(
sm = [0.11, 0.12, 0.13, 0.09],
d = [ 5, 10, 20, 40],
bd = [ 1.5, 1.5, 1.5, 1.5],
r = 55
),
Profile(
sm = [0.31, 0.33, 0.33, 0.49],
d = [ 5, 10, 20, 40],
bd = [ 1.5, 1.5, 1.5, 1.5],
r = 95
)
]
print("Set up %.0f profiles." % len(profiles))
# %%
i = 0
for P in profiles:
i += 1
# First order estimate of the average soil moisture
sm_estimate = P.sm.mean()
# Neutron penetration depth at this location
P.D86 = D86(
sm = P.sm.mean(),
bd = P.bd.mean(),
r = P.r
)
# Weights for this profile
# data = P.data.copy()
P.data["weight"] = Weight_d(
P.data["d"],
P.D86
)
# Calculate weighted sm average
P.sm_wavg = np.average(
P.data["sm"],
weights = P.data["weight"]
)
# Calculate horizontal weight for this profile
P.w_r = Wr(
r = P.r,
sm = P.sm_wavg,
hum= air_humidity
)
print("Profile %.0f: vertical wt. avg. SM: %.3f, CRNS depth: %3.0f cm, distance: %2.0f m, horizontal weight: %7.0f" \
% (i, P.sm_wavg, P.D86, P.r, P.w_r))
# %%
# Horizontal average
profiles_sm_wavg = [P.sm_wavg for P in profiles]
profiles_weights = [P.w_r for P in profiles]
sm_wavg = np.average(
profiles_sm_wavg,
weights = profiles_weights
)
print("Horizontal wt. avg. SM: %.3f" % sm_wavg)
# %%
# Footprint
footprint_m = get_footprint(sm_wavg, air_humidity, air_pressure)
print("Footprint radius: %.0f m" % footprint_m)
# %%
# Plot vertical profiles
plt.title("SM profiles, vertical wt. avg. SM, and CRNS depth")
for P in profiles:
plt.plot(P.sm, P.d, "o-")
plt.plot(
[P.sm_wavg, P.sm_wavg],
[0, P.D86],
"v--", color="black", ms=10
)
plt.xlim(0,1)
plt.ylim(60,0)
plt.xlabel("SM")
plt.ylabel("Depth (cm)")
plt.show()
# %%
# Plot horizontal profiles
plt.title("SM profile locations and vertical wt. avg. SM,\nhorizontal wt. avg. SM, CRNS footprint")
for P in profiles:
plt.plot(P.r, P.sm_wavg, "o-")
plt.plot(
[0, footprint_m],
[sm_wavg, sm_wavg],
">--", color="black", ms=10
)
plt.xlim(0,250)
plt.ylim(0.0,1.0)
plt.xlabel("Distance (m)")
plt.ylabel("SM")
plt.show()
# %%
examples/weighting/Schroen2017hess.py 0 → 100644
+ 95
− 0
View file @ ab62011a
"""
CoRNy CORN
Functions specific for data processing of COsmic-Ray Neutrons
from: https://git.ufz.de/CRNS/cornish_pasdy/-/blob/master/corny/corn.py
version: 0.62
"""
import pandas
import pandas as pd
import numpy as np
import os
def app_file_path(path):
return(path)
thisfolder = os.path.dirname(__file__)
# Footprint
def D86(sm, bd=1, r=1):
return(1/bd*( 8.321+0.14249*( 0.96655+np.exp(-r/100))*(20.0+sm) / (0.0429+sm)))
def Weight_d(d, D):
return(np.exp(-2*d/D))
# Preload footprint radius matrix
preload_footprint_data = np.loadtxt(app_file_path('footprint_radius.csv'))
# preload_footprint_data = np.loadtxt(os.path.join(thisfolder, 'footprint_radius.csv'))
def get_footprint(sm, h, p, lookup_file=None):
"""
usage: data['footprint_radius'] = data.apply(lambda row: get_footprint( row[smv], row[ah], row[p] ), axis=1)
"""
if np.isnan(sm) or np.isnan(h) or np.isnan(p): return(np.nan)
if sm < 0.01: sm = 0.01
if sm > 0.49: sm = 0.49
if h > 30 : h = 30
#print(sm, h, p, int(round(100*sm)), int(round(h)))
if lookup_file is None:
footprint_data = preload_footprint_data
#footprint_data = np.loadtxt(os.path.join(pasdy_path ,'corny/footprint_radius.csv'))
elif isinstance(lookup_file, str):
if os.path.exists(lookup_file):
footprint_data = np.loadtxt(lookup_file)
else:
return(np.nan)
elif isinstance(lookup_file, np.ndarray):
footprint_data = lookup_file
else:
return(np.nan)
return(footprint_data[int(round(100*sm))][int(round(h))] * 0.4922/(0.86-np.exp(-p/1013.25)))
def get_footprint_volume(depth, radius, theta, bd):
return((depth + D86(theta, bd, radius))*0.01*0.47*radius**2*3.141 /1000)
# 0.44 (dry) ..0.5 (wet) is roughly the average D over radii
def Wr_approx(r=1):
return((30*np.exp(-r/1.6)+np.exp(-r/100))*(1-np.exp(-3.7*r)))
def Wr(r=1, sm=0.1, hum=5, normalize=False):
x = hum
y = sm
a00 = 8735; a01 = 22.689; a02 = 11720; a03 = 0.00978; a04 = 9306; a05 = 0.003632
a10 = 2.7925e-002; a11 = 6.6577; a12 = 0.028544; a13 = 0.002455; a14 = 6.851e-005; a15 = 12.2755
a20 = 247970; a21 = 23.289; a22 = 374655; a23 = 0.00191; a24 = 258552
a30 = 5.4818e-002; a31 = 21.032; a32 = 0.6373; a33 = 0.0791; a34 = 5.425e-004
b00 = 39006; b01 = 15002337; b02 = 2009.24; b03 = 0.01181; b04 = 3.146; b05 = 16.7417; b06 = 3727
b10 = 6.031e-005; b11 = 98.5; b12 = 0.0013826
b20 = 11747; b21 = 55.033; b22 = 4521; b23 = 0.01998; b24 = 0.00604; b25 = 3347.4; b26 = 0.00475
b30 = 1.543e-002; b31 = 13.29; b32 = 1.807e-002; b33 = 0.0011; b34 = 8.81e-005; b35 = 0.0405; b36 = 26.74
A0 = (a00*(1+a03*x)*np.exp(-a01*y)+a02*(1+a05*x)-a04*y)
A1 = ((-a10+a14*x)*np.exp(-a11*y/(1+a15*y))+a12)*(1+x*a13)
A2 = (a20*(1+a23*x)*np.exp(-a21*y)+a22-a24*y)
A3 = a30*np.exp(-a31*y)+a32-a33*y+a34*x
B0 = (b00-b01/(b02*y+x-0.13))*(b03-y)*np.exp(-b04*y)-b05*x*y+b06
B1 = b10*(x+b11)+b12*y
B2 = (b20*(1-b26*x)*np.exp(-b21*y*(1-x*b24))+b22-b25*y)*(2+x*b23)
B3 = ((-b30+b34*x)*np.exp(-b31*y/(1+b35*x+b36*y))+b32)*(2+x*b33)
if np.isscalar(r):
if r <= 1: w = (A0*(np.exp(-A1*r)) + A2*np.exp(-A3*r))*(1-np.exp(-3.7*r))
elif (r > 1) & (r < 50): w = A0*(np.exp(-A1*r)) + A2*np.exp(-A3*r)
elif (r >= 50): w = B0*(np.exp(-B1*r)) + B2*np.exp(-B3*r)
return(w)
else:
W = pandas.DataFrame()
W['r'] = r
W['w'] = 0
W.loc[W.r <= 1,'w'] = (A0*(np.exp(-A1*W.loc[W.r <= 1,'r'])) + A2*np.exp(-A3*W.loc[W.r <= 1,'r']))*(1-np.exp(-3.7*W.loc[W.r <= 1,'r']))
W.loc[W.r > 1,'w'] = A0*(np.exp(-A1*W.loc[W.r > 1,'r'])) + A2*np.exp(-A3*W.loc[W.r > 1,'r'])
W.loc[W.r >= 50,'w'] = B0*(np.exp(-B1*W.loc[W.r >= 50,'r'])) + B2*np.exp(-B3*W.loc[W.r >= 50,'r'])
if normalize:
W.w /= W.w.sum()
return(W.w.values)
examples/weighting/footprint_radius.csv 0 → 100644
+ 50
− 0
View file @ ab62011a
231 229 227 225 223 221 219 217 214 212 210 207 205 202 200 197 195 193 190 188 186 183 181 179 177 175 172 170 168 166
226 224 222 219 217 214 212 210 207 205 203 201 198 196 194 192 190 188 186 184 181 179 177 175 174 172 170 168 166 164
229 226 223 220 218 215 212 210 208 205 203 200 198 196 194 192 189 187 185 183 181 179 177 175 173 171 170 168 166 164
231 228 225 222 220 217 214 212 209 206 204 202 199 197 195 192 190 188 186 184 182 180 178 176 174 172 170 168 167 165
233 230 227 224 221 218 215 213 210 207 205 203 200 198 195 193 191 189 187 185 182 180 178 177 175 173 171 169 167 165
233 230 227 224 221 218 216 213 210 208 205 203 200 198 196 193 191 189 187 185 183 181 179 177 175 173 171 169 167 166
233 229 226 223 221 218 215 212 210 207 205 202 200 197 195 193 191 188 186 184 182 180 178 176 174 172 171 169 167 165
231 228 225 222 219 217 214 211 209 206 203 201 199 196 194 192 190 187 185 183 181 179 177 175 173 172 170 168 166 164
229 226 223 220 217 215 212 209 207 204 202 199 197 195 192 190 188 186 184 182 180 178 176 174 172 170 168 167 165 163
226 223 221 218 215 212 210 207 204 202 200 197 195 193 190 188 186 184 182 180 178 176 174 172 170 168 167 165 163 161
223 220 218 215 212 209 207 204 202 199 197 195 192 190 188 186 184 181 179 177 175 174 172 170 168 166 164 163 161 159
220 217 214 212 209 206 204 201 199 196 194 192 190 187 185 183 181 179 177 175 173 171 169 167 166 164 162 160 159 157
217 214 211 208 206 203 201 198 196 193 191 189 186 184 182 180 178 176 174 172 170 168 167 165 163 161 160 158 156 155
213 210 207 205 202 200 197 195 192 190 188 185 183 181 179 177 175 173 171 169 167 166 164 162 160 159 157 155 154 152
209 206 204 201 198 196 194 191 189 187 184 182 180 178 176 174 172 170 168 166 164 163 161 159 157 156 154 153 151 149
205 202 200 197 195 192 190 188 185 183 181 179 177 175 173 171 169 167 165 163 161 160 158 156 155 153 151 150 148 147
201 198 196 193 191 189 186 184 182 179 177 175 173 171 169 167 165 164 162 160 158 156 155 153 152 150 148 147 145 144
197 194 192 189 187 185 182 180 178 176 174 172 170 168 166 164 162 160 158 157 155 153 152 150 149 147 145 144 142 141
193 190 188 186 183 181 179 176 174 172 170 168 166 164 162 161 159 157 155 154 152 150 149 147 146 144 143 141 140 138
189 186 184 182 179 177 175 173 171 169 167 165 163 161 159 157 155 154 152 150 149 147 146 144 143 141 140 138 137 135
185 182 180 178 176 173 171 169 167 165 163 161 159 158 156 154 152 151 149 147 146 144 143 141 140 138 137 135 134 133
181 179 176 174 172 170 168 166 164 162 160 158 156 154 152 151 149 147 146 144 143 141 140 138 137 135 134 133 131 130
177 175 172 170 168 166 164 162 160 158 156 154 153 151 149 148 146 144 143 141 140 138 137 135 134 132 131 130 129 127
173 171 169 167 165 163 161 159 157 155 153 151 149 148 146 144 143 141 140 138 137 135 134 132 131 130 128 127 126 125
170 167 165 163 161 159 157 155 153 152 150 148 146 145 143 141 140 138 137 135 134 133 131 130 128 127 126 125 123 122
166 164 162 160 158 156 154 152 150 148 147 145 143 142 140 139 137 136 134 133 131 130 128 127 126 125 123 122 121 120
163 161 158 156 155 153 151 149 147 145 144 142 140 139 137 136 134 133 131 130 129 127 126 125 123 122 121 120 119 117
159 157 155 153 151 150 148 146 144 143 141 139 138 136 135 133 132 130 129 127 126 125 124 122 121 120 119 117 116 115
156 154 152 150 148 147 145 143 141 140 138 137 135 134 132 131 129 128 126 125 124 122 121 120 119 118 116 115 114 113
153 151 149 148 146 144 142 140 139 137 136 134 133 131 130 128 127 125 124 123 121 120 119 118 117 115 114 113 112 111
150 149 147 145 143 141 140 138 136 135 133 132 130 129 127 126 125 123 122 121 119 118 117 116 115 113 112 111 110 109
148 146 144 142 141 139 137 136 134 132 131 129 128 127 125 124 122 121 120 119 117 116 115 114 113 111 110 109 108 107
145 143 142 140 138 137 135 133 132 130 129 127 126 124 123 122 120 119 118 117 115 114 113 112 111 110 109 108 107 105
143 141 139 138 136 134 133 131 130 128 127 125 124 123 121 120 119 117 116 115 114 112 111 110 109 108 107 106 105 104
141 139 137 136 134 132 131 129 128 126 125 123 122 121 119 118 117 116 114 113 112 111 110 109 108 106 105 104 103 102
139 137 135 134 132 131 129 127 126 125 123 122 120 119 118 116 115 114 113 112 110 109 108 107 106 105 104 103 102 101
137 135 134 132 130 129 127 126 124 123 122 120 119 118 116 115 114 113 111 110 109 108 107 106 105 104 103 102 101 100
136 134 132 131 129 127 126 124 123 122 120 119 117 116 115 114 112 111 110 109 108 107 106 105 104 102 101 100 100 99
134 132 131 129 128 126 124 123 122 120 119 117 116 115 114 112 111 110 109 108 107 106 104 103 102 101 100 99 98 97
133 131 129 128 126 125 123 122 120 119 118 116 115 114 113 111 110 109 108 107 106 104 103 102 101 100 99 98 97 97
132 130 128 127 125 124 122 121 119 118 117 115 114 113 112 110 109 108 107 106 105 104 102 101 100 99 98 98 97 96
131 129 127 126 124 123 121 120 118 117 116 114 113 112 111 109 108 107 106 105 104 103 102 101 100 99 98 97 96 95
130 128 127 125 123 122 120 119 118 116 115 114 112 111 110 109 108 106 105 104 103 102 101 100 99 98 97 96 95 94
129 127 126 124 123 121 120 118 117 116 114 113 112 110 109 108 107 106 105 104 102 101 100 99 98 97 96 95 94 94
129 127 125 124 122 121 119 118 116 115 114 112 111 110 109 108 106 105 104 103 102 101 100 99 98 97 96 95 94 93
128 127 125 123 122 120 119 117 116 115 113 112 111 110 108 107 106 105 104 103 101 100 99 98 97 96 95 94 94 93
128 126 125 123 122 120 119 117 116 114 113 112 110 109 108 107 106 104 103 102 101 100 99 98 97 96 95 94 93 92
128 126 124 123 121 120 118 117 116 114 113 112 110 109 108 107 105 104 103 102 101 100 99 98 97 96 95 94 93 92
128 126 124 123 121 120 118 117 115 114 113 111 110 109 108 106 105 104 103 102 101 100 99 98 97 96 95 94 93 92
128 126 124 123 121 120 118 117 115 114 113 111 110 109 108 106 105 104 103 102 101 100 99 97 96 95 94 94 93 92
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment