Commit 57f48da2 authored by Niklas Bohn's avatar Niklas Bohn

Tested different version of new surface model.

parent 4efa0d09
......@@ -128,7 +128,8 @@ class FoGen(object):
self.pt[:, :, 3] = np.full(self.data.shape[:2], options["metadata"]["aot"])
self.pt[:, :, 4] = self.raa
self.par_opt = ["cwv", "cwc", "ice", "lai", "dasf", "p_max", "k", "b"]
# self.par_opt = ["cwv", "cwc", "ice", "lai", "dasf", "p_max", "k", "b"]
self.par_opt = ["cwv", "ewt", "ice", "dasf", "p"]
self.state_dim = len(self.par_opt)
self.prior_mean, self.use_prior_mean, self.ll, self.ul, self.prior_sigma, self.unknowns = [], [], [], [], [], []
......@@ -257,13 +258,13 @@ class FoGen(object):
w0 = np.exp(-xx[1] * 1e7 * self.abs_co_w - xx[2] * 1e7 * self.abs_co_i)
# canopy photon re-collision probability (Smolander & Stenberg 2005)
p = xx[5] * (1 - np.exp(-xx[6] * xx[3] ** xx[7]))
# p = xx[5] * (1 - np.exp(-xx[6] * xx[3] ** xx[7]))
# canopy single scattering albedo (Lewis & Disney 2007)
W = w0 * (1 - p) / (1 - p * w0)
W = w0 * (1 - xx[4]) / (1 - xx[4] * w0)
# canopy bidirectional reflectance factor (Knyazikhin et al. 2013)
rho = xx[4] * W
rho = xx[3] * W
return rho
......@@ -899,7 +900,7 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
# perform first guess liquid water retrieval based on the NDWI
if fo.use_prior_mean[1]:
cwc_fg = np.full(fo.data.shape[:-2], fo.prior_mean[1])
ewt_fg = np.full(fo.data.shape[:2], fo.prior_mean[1])
else:
logger.info("perform first guess liquid water retrieval based on the NDWI...")
ndwi = np.zeros(fo.data.shape[:2])
......@@ -907,7 +908,7 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
for jj in range(fo.data.shape[1]):
ndwi[ii, jj] = (fo.data[ii, jj, 72] - fo.data[ii, jj, 118]) / \
(fo.data[ii, jj, 72] + fo.data[ii, jj, 118])
cwc_fg = 0.49585423 * ndwi**2 - 0.0396154 * ndwi - 0.0174602
ewt_fg = 0.49585423 * ndwi**2 - 0.0396154 * ndwi - 0.0174602
# perform first guess snow retrieval based on the NDSI
if fo.use_prior_mean[2]:
......@@ -936,9 +937,9 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
for ii in range(xa.shape[0]):
for jj in range(xa.shape[1]):
xa[ii, jj, 0] = cwv_fg[ii, jj]
xa[ii, jj, 1] = cwc_fg[ii, jj]
xa[ii, jj, 1] = ewt_fg[ii, jj]
xa[ii, jj, 2] = ice_fg[ii, jj]
for dd in range(fo.state_dim - 5, fo.state_dim):
for dd in range(fo.state_dim - 2, fo.state_dim):
xa[ii, jj, dd] = fo.prior_mean[dd]
for dd in range(fo.state_dim):
......@@ -965,13 +966,13 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
globs["__forward__"] = opt_func
globs["__cwv_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__cwc_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__ewt_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__ice_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__lai_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
# globs["__lai_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__dasf_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__p_max_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__k_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__b_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__p_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
# globs["__k_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
# globs["__b_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]))
globs["__toa_model__"] = SharedNdarray(dims=list(fo.data.shape[:2]) + [fo.fit_wvl.shape[0]])
globs["__sx__"] = SharedNdarray(dims=list(fo.data.shape[:2]) + [fo.state_dim] + [fo.state_dim])
......@@ -1013,13 +1014,13 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
t1 = time()
cwv_model = globs["__cwv_model__"].np
cwc_model = globs["__cwc_model__"].np
ewt_model = globs["__ewt_model__"].np
ice_model = globs["__ice_model__"].np
lai_model = globs["__lai_model__"].np
# lai_model = globs["__lai_model__"].np
dasf_model = globs["__dasf_model__"].np
p_max_model = globs["__p_max_model__"].np
k_model = globs["__k_model__"].np
b_model = globs["__b_model__"].np
p_model = globs["__p_model__"].np
# k_model = globs["__k_model__"].np
# b_model = globs["__b_model__"].np
toa_model = globs["__toa_model__"].np
sx = globs["__sx__"].np
......@@ -1029,9 +1030,8 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
logger.info("Done!")
logger.info("Runtime: %.2f" % (t1 - t0) + " s")
res = {"cwv_model": cwv_model, "cwc_model": cwc_model, "ice_model": ice_model, "lai_model": lai_model,
"dasf_model": dasf_model, "p_max_model": p_max_model, "k_model": k_model, "b_model": b_model,
"toa_model": toa_model, "sx": sx, "scem": scem, "srem": srem}
res = {"cwv_model": cwv_model, "ewt_model": ewt_model, "ice_model": ice_model, "dasf_model": dasf_model,
"p_model": p_model, "toa_model": toa_model, "sx": sx, "scem": scem, "srem": srem}
return res
......@@ -1053,13 +1053,13 @@ def __oe__(ii):
model = __forward__(xx=res[0], pt=__pt__[i1, i2, :], dt=__data__[i1, i2, :], model_output=True)
__cwv_model__[i1, i2] = res[0][0]
__cwc_model__[i1, i2] = res[0][1]
__ewt_model__[i1, i2] = res[0][1]
__ice_model__[i1, i2] = res[0][2]
__lai_model__[i1, i2] = res[0][__state_dim__ - 5]
__dasf_model__[i1, i2] = res[0][__state_dim__ - 4]
__p_max_model__[i1, i2] = res[0][__state_dim__ - 3]
__k_model__[i1, i2] = res[0][__state_dim__ - 2]
__b_model__[i1, i2] = res[0][__state_dim__ - 1]
# __lai_model__[i1, i2] = res[0][__state_dim__ - 5]
__dasf_model__[i1, i2] = res[0][__state_dim__ - 2]
__p_model__[i1, i2] = res[0][__state_dim__ - 1]
# __k_model__[i1, i2] = res[0][__state_dim__ - 2]
# __b_model__[i1, i2] = res[0][__state_dim__ - 1]
__toa_model__[i1, i2, :] = model
# a posteriori covariance matrix
......
......@@ -30,56 +30,35 @@
"use_prior_mean": false,
"ll": 0.0001,
"ul": 4.9999,
"prior_sigma": 1000
"prior_sigma": 0.5
},
"cwc": {
"prior_mean": 0.02,
"use_prior_mean": false,
"ewt": {
"prior_mean": 0.0111,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.4999,
"prior_sigma": 0.01
"ul": 0.0999,
"prior_sigma": 0.006979349
},
"ice": {
"prior_mean": 0,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.4999,
"prior_sigma": 3.1622776601683795e-05
},
"lai": {
"prior_mean": 5,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 12.9999,
"prior_sigma": 1000
"prior_sigma": 1e-100
},
"dasf": {
"prior_mean": 0.25,
"prior_mean": 0.15,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.9999,
"prior_sigma": 1000
"prior_sigma": 0.000001
},
"p_max": {
"prior_mean": 0.88,
"p": {
"prior_mean": 0.7954979526021773,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.9999,
"prior_sigma": 1000
},
"k": {
"prior_mean": 0.7,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.9999,
"prior_sigma": 1000
},
"b": {
"prior_mean": 0.75,
"use_prior_mean": true,
"ll": 0.0001,
"ll": 0.4001,
"ul": 0.9999,
"prior_sigma": 1000
"prior_sigma": 0.13466832630044967
}
},
"unknowns": {
......
{
"sensor":{
"name": "EnMAP",
"fit":{
"wvl_center": [1051.3, 1062.6, 1074, 1085.4, 1096.9, 1108.5, 1120.1, 1131.8, 1143.5, 1155.3, 1167.1, 1179,
1190.9, 1202.8, 1214.8, 1226.7, 1238.7, 1250.7],
"fwhm": [13.52, 13.62, 13.72, 13.8, 13.88, 13.96, 14.03, 14.09, 14.15, 14.19, 14.24, 14.28, 14.31, 14.34,
14.36, 14.38, 14.39, 14.4],
"snr": [322.94254006, 313.45752342, 313.32464722, 306.19455173, 281.35641378, 227.9759974, 157.22390595,
148.93620623, 154.40197388, 182.60431866, 232.89644996, 250.1244036, 252.2267179, 272.16592448,
299.71964816, 316.11909184, 326.33202741, 312.28461288],
"idx": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
}
},
"metadata":{
"jday": 16,
"month": 6,
"sza": 40,
"saa": 360,
"vza": 0,
"hsf": 0,
"aot": 0.2
},
"retrieval": {
"fn_LUT": "",
"cpu":12,
"state_vector": {
"cwv": {
"prior_mean": 2.5,
"use_prior_mean": false,
"ll": 0.0001,
"ul": 4.9999,
"prior_sigma": 1000
},
"ewt": {
"prior_mean": 0.02,
"use_prior_mean": false,
"ll": 0.0001,
"ul": 0.4999,
"prior_sigma": 0.01
},
"ice": {
"prior_mean": 0,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.4999,
"prior_sigma": 3.1622776601683795e-05
},
"lai": {
"prior_mean": 5,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 12.9999,
"prior_sigma": 1000
},
"dasf": {
"prior_mean": 0.25,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.9999,
"prior_sigma": 1000
},
"p": {
"prior_mean": 0.7,
"use_prior_mean": true,
"ll": 0.4001,
"ul": 0.9999,
"prior_sigma": 1000
},
"p_max": {
"prior_mean": 0.88,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.9999,
"prior_sigma": 1000
},
"k": {
"prior_mean": 0.7,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.9999,
"prior_sigma": 1000
},
"b": {
"prior_mean": 0.75,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.9999,
"prior_sigma": 1000
}
},
"unknowns": {
"Skyview": {
"sigma": 0.1
},
"H2O_ABSCO": {
"sigma": 0.01
},
"Liquid_ABSCO": {
"sigma": 2.43e-07
},
"Ice_ABSCO": {
"sigma": 2.39e-07
}
},
"inversion": {
"gnform": "n",
"full": false,
"maxiter": 35,
"eps": 0.01
}
}
}
......@@ -55,10 +55,9 @@ def make_ac_generic(data_l1b, fo, xx, logger=None):
for irow in range(nrows):
data_l2a[irow, icol, :] = fo.surf_ref(dt=data_l1b[irow, icol, :], xx=[xx["cwv_model"][irow, icol]],
pt=fo.pt[irow, icol, :], mode="full")
swir_feature = fo.surface_model(xx=[xx["cwv_model"][irow, icol], xx["cwc_model"][irow, icol],
xx["ice_model"][irow, icol], xx["lai_model"][irow, icol],
xx["dasf_model"][irow, icol], xx["p_max_model"][irow, icol],
xx["k_model"][irow, icol], xx["b_model"][irow, icol]])
swir_feature = fo.surface_model(xx=[xx["cwv_model"][irow, icol], xx["ewt_model"][irow, icol],
xx["ice_model"][irow, icol], xx["dasf_model"][irow, icol],
xx["p_model"][irow, icol]])
data_l2a[irow, icol, fo.fit_wvl] = swir_feature
return data_l2a
......
Markdown is supported
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