Skip to content

Commit

Permalink
add conversion from genuchten to brooks
Browse files Browse the repository at this point in the history
  • Loading branch information
@martinvonk
martinvonk committed Dec 5, 2024
1 parent b6cf6d3 commit a15a1a2
Showing 1 changed file with 158 additions and 38 deletions.
196 changes: 158 additions & 38 deletions src/pedon/soildata.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@
# "pandas",
# "geopandas",
# "matplotlib",
# "openpyxl"
# "openpyxl",
# "pedon",
# ]
# ///
import requests
Expand All @@ -14,60 +15,72 @@
from shapely import Polygon


#%%
# %%
# voor 1 punt
# send get request
r = requests.get(
url = "https://www.soilphysics.wur.nl/soil.php",
params = {"latitude": 52, "longitude": 5.5} )
url="https://www.soilphysics.wur.nl/soil.php",
params={"latitude": 52, "longitude": 5.5},
)

# extract data in json format
data = r.json()

# further processing
pd.DataFrame(data["horizon"])

#%%
# %%
# voor bofek
# haal cluster en bijbehorende profielen op
data = pd.read_excel('tabel_BOFEK_2020_V_1_0.xlsx', sheet_name="Data (2)", skiprows=9, index_col=0)
profiel_cluster = dict(zip(data['Cluster'], data["Profiel"]))
#%%
data = pd.read_excel(
"tabel_BOFEK_2020_V_1_0.xlsx", sheet_name="Data (2)", skiprows=9, index_col=0
)
profiel_cluster = dict(zip(data["Cluster"], data["Profiel"]))
# %%
# haal de geometrie van de clusters op
gdf = gpd.read_file('BOFEK2020.gdb')
gdf = gpd.read_file("BOFEK2020.gdb")
# snij bij met extent van wijde wormer
extent = Polygon([(117500.0, 497250.0), (125000.0, 497250.0), (125000.0, 503000.0), (117500.0, 503000.0)]) # wijde wormer extent
gdf_clip = gpd.clip(gdf, extent).dropna(subset=['BOFEK2020'])
extent = Polygon(
[
(117500.0, 497250.0),
(125000.0, 497250.0),
(125000.0, 503000.0),
(117500.0, 503000.0),
]
) # wijde wormer extent
gdf_clip = gpd.clip(gdf, extent).dropna(subset=["BOFEK2020"])

#%%
# %%
# pak unieke profielen op basis van clusters
unique_profiles = gdf_clip['BOFEK2020'].astype(int).map(profiel_cluster).unique()
#%%
unique_profiles = gdf_clip["BOFEK2020"].astype(int).map(profiel_cluster).unique()
# %%
# haal de profielen op
profiles = pd.read_csv(
"AllProfiles_368.csv",
index_col=0,
na_values=99999,
dtype = {'iSoil1': 'int64',
'iSoil2': 'int64',
'iSoil3': 'int64',
'iSoil4': 'int64',
'iSoil5': 'int64',
'iSoil6': 'int64',
'iSoil7': 'int64',
'iSoil8': 'int64',
'iSoil9': 'int64',
'iZ1': 'float64',
'iZ2': 'float64',
'iZ3': 'float64',
'iZ4': 'float64',
'iZ5': 'float64',
'iZ6': 'float64',
'iZ7': 'float64',
'iZ8': 'float64',
'iZ9': 'float64'},
).drop(columns=["iBOFEK2012", "iHoofd"])
#%%
dtype={
"iSoil1": "int64",
"iSoil2": "int64",
"iSoil3": "int64",
"iSoil4": "int64",
"iSoil5": "int64",
"iSoil6": "int64",
"iSoil7": "int64",
"iSoil8": "int64",
"iSoil9": "int64",
"iZ1": "float64",
"iZ2": "float64",
"iZ3": "float64",
"iZ4": "float64",
"iZ5": "float64",
"iZ6": "float64",
"iZ7": "float64",
"iZ8": "float64",
"iZ9": "float64",
},
).drop(columns=["iBOFEK2012", "iHoofd"])
# %%
# haal de bodemparameters op (kan netter)
soil_mapper = {
1: pe.Soil("B01").from_staring("2018").model,
Expand Down Expand Up @@ -108,15 +121,122 @@
36: pe.Soil("O18").from_staring("2018").model,
}

#%%
# zet bofek profiel om te zetten naar boring met bodemparameters
# %%
# zet Bofek profiel om te zetten naar boring met bodemparameters


def bofek_to_boring(profiles: pd.DataFrame, iprofile: int):
profile = profiles.loc[[iprofile]]
depth = profile.loc[:, ["iZ1", "iZ2", "iZ3", "iZ4", "iZ5", "iZ6", "iZ7", "iZ8", "iZ9"]].squeeze().rename("Depth [cm]")
soil = profile.loc[:, ["iSoil1", "iSoil2", "iSoil3", "iSoil4", "iSoil5", "iSoil6", "iSoil7", "iSoil8", "iSoil9"]].squeeze().rename("Soil")
depth = (
profile.loc[:, ["iZ1", "iZ2", "iZ3", "iZ4", "iZ5", "iZ6", "iZ7", "iZ8", "iZ9"]]
.squeeze()
.rename("Depth [cm]")
)
soil = (
profile.loc[
:,
[
"iSoil1",
"iSoil2",
"iSoil3",
"iSoil4",
"iSoil5",
"iSoil6",
"iSoil7",
"iSoil8",
"iSoil9",
],
]
.squeeze()
.rename("Soil")
)
soil.index = depth
# soil = soil[~soil.duplicated(keep="last")] # remove duplicates maar dit klopt zo nog niet als je meerdere trajecten hebt in bodem met dezelfde bodem
return soil.map(soil_mapper).dropna().rename(iprofile)


{iprofile: bofek_to_boring(profiles, iprofile) for iprofile in unique_profiles}
# %%
import pedon as pe

import numpy as np


gs = [
pe.Genuchten(0.1, 0.05, 0.4, 0.3, 1.2),
pe.Genuchten(0.1, 0.05, 0.4, 0.3, 1.2),
pe.Genuchten(0.1, 0.05, 0.4, 0.3, 1.2),
]

depth = np.array([0.0, 20.0, 40.0, 70.0, 120.0])
ghg_depth = 50.0
h = depth - ghg_depth # pressure head
# pressure head = head - z = h + depths
dh = 1.0
for hs, he, g in zip(h[:-1], h[1:], gs):
hr = np.arange(hs, he, dh)
print(hr)
theta = g.theta(hr)


# %%
def get_brooks(gen: pe.Genuchten, h: np.ndarray[float] | None) -> pe.Brooks:
if h is not None:
k = gen.k(h)
theta = gen.theta(h)
soilsample = pe.SoilSample(h=h, k=k, theta=theta)
pbounds = pe._params.pBrooks.copy()
pbounds.loc["theta_r"] = (
gen.theta_r,
max(g.theta_r - 0.1, 0.0),
gen.theta_r + 0.1,
)
pbounds.loc["theta_s"] = (gen.theta_s, gen.theta_s - 0.1, gen.theta_s + 0.1)
pbounds.loc["l"] = (5.0, 0.01, 20.0)
pbounds.loc["h_b"] = (1.0, 0.1, 100.0)
bc = soilsample.fit(pe.Brooks, pbounds=pbounds, k_s=g.k_s)
else:
# Morel-Seytoux (1996) - Parameter equivalence for the Brooks-Corey and van Genuchten soil characteristics
eps = 1 + 2 / gen.m # eq 16b
h_b = (
1
/ gen.alpha
* (eps + 3)
/ (2 * eps * (eps - 1))
* (147.8 + 8.1 * eps + 0.0928 * eps**2)
/ (55.6 + 7.4 * eps + eps**2)
) # eq 17
l = 2 / (eps - 3) # because eps = 3 + 2 / l
bc = pe.Brooks(
k_s=gen.k_s, theta_r=gen.theta_r, theta_s=gen.theta_s, h_b=h_b, l=l
)
return bc


# %%

# name = "B05"
# g = pe.Soil(name=name).from_staring("2018").model
hb = []
for name, g in soil_mapper.items():
h = np.logspace(-2, 6, num=11)
bc1 = get_brooks(g, h)
bc2 = get_brooks(g, None)

print(name)
print(3 + 2 / bc1.l, 3 + 2 / bc2.l)
print(bc1.h_b, bc2.h_b)

ax1 = pe.plot_swrc(g)
pe.plot_swrc(bc1, ax=ax1)
pe.plot_swrc(bc2, ax=ax1)
ax1.set_title(name)

ax2 = pe.plot_hcf(g)
pe.plot_hcf(bc1, ax=ax2)
pe.plot_hcf(bc2, ax=ax2)
ax2.set_title(name)
hb.append(bc1.h_b)
# break

# %%

0 comments on commit a15a1a2

Please sign in to comment.