Skip to content

API Reference

Core Models

Soil profile model representing a vertical arrangement of soil layers.

SoilProfile

Bases: BaseModel

A soil profile composed of layers with spatial information.

Source code in simplesoilprofile/models/profile.py 
 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
class SoilProfile(BaseModel):
    """A soil profile composed of layers with spatial information."""

    model_config = {
        'arbitrary_types_allowed': True,
    }

    name: str = Field(..., description="Name or identifier of the soil profile")
    description: str | None = Field(None, description="Description of the soil profile")

    location: Point | None = Field(None, description="Spatial location of the profile as a Point object (x, y coordinates)")
    elevation: float | None = Field(None, description="Z coordinate (elevation) of the profile surface [m]")

    layers: list[SoilLayer] = Field(..., description="List of soil layers in the profile")
    layer_bottoms: list[float] = Field(
        ...,
        description="List of bottom depths for each layer [cm]. Top of first layer is assumed to be 0."
    )

    @model_validator(mode='after')
    def validate_layer_depths(self) -> 'SoilProfile':
        """Validate that layer depths are consistent and monotonically increasing."""
        if not self.layer_bottoms:
            raise ValueError("Layer depths cannot be empty")

        if len(self.layer_bottoms) != len(self.layers):
            raise ValueError("Must provide bottom depths for all layers")

        prev_depth = 0
        for i, bottom in enumerate(self.layer_bottoms):
            if bottom <= prev_depth:
                raise ValueError(f"Layer {i}: bottom depth must be greater than previous layer")
            prev_depth = bottom

        return self

    @computed_field
    @property
    def layer_bounds(self) -> dict[int, tuple[float, float]]:
        """Get dictionary of (top, bottom) bounds for each layer."""
        bounds = {}
        prev_depth = 0
        for i, bottom in enumerate(self.layer_bottoms):
            bounds[i] = (prev_depth, bottom)
            prev_depth = bottom
        return bounds

    @computed_field
    @property
    def profile_depth(self) -> float:
        """Total depth of the profile in cm (bottom of last layer)."""
        return self.layer_bottoms[-1] if self.layer_bottoms else 0.0

    def get_layer_at_depth(self, depth: float) -> tuple[SoilLayer, int] | None:
        """Get the soil layer at a specific depth."""
        prev_depth = 0
        for i, bottom in enumerate(self.layer_bottoms):
            if prev_depth <= depth <= bottom:
                return (self.layers[i], i)
            prev_depth = bottom
        return None

    def get_sublayer_boundaries(self) -> dict[int, list[float]]:
        """Get all sublayer boundaries for each layer in the profile."""
        boundaries = {}
        prev_depth = 0
        for i, bottom in enumerate(self.layer_bottoms):
            boundaries[i] = self.layers[i].get_sublayer_boundaries(prev_depth, bottom)
            prev_depth = bottom
        return boundaries

    def get_sublayer_depths(self) -> list[float]:
        """Get a sorted list of all sublayer boundary depths in the profile."""
        all_depths = []
        for depths in self.get_sublayer_boundaries().values():
            all_depths.extend(depths)
        return sorted(set(all_depths))

layer_bounds property

Get dictionary of (top, bottom) bounds for each layer.

profile_depth property

Total depth of the profile in cm (bottom of last layer).

get_layer_at_depth(depth)

Get the soil layer at a specific depth.

Source code in simplesoilprofile/models/profile.py 
62
63
64
65
66
67
68
69
def get_layer_at_depth(self, depth: float) -> tuple[SoilLayer, int] | None:
    """Get the soil layer at a specific depth."""
    prev_depth = 0
    for i, bottom in enumerate(self.layer_bottoms):
        if prev_depth <= depth <= bottom:
            return (self.layers[i], i)
        prev_depth = bottom
    return None

get_sublayer_boundaries()

Get all sublayer boundaries for each layer in the profile.

Source code in simplesoilprofile/models/profile.py 
71
72
73
74
75
76
77
78
def get_sublayer_boundaries(self) -> dict[int, list[float]]:
    """Get all sublayer boundaries for each layer in the profile."""
    boundaries = {}
    prev_depth = 0
    for i, bottom in enumerate(self.layer_bottoms):
        boundaries[i] = self.layers[i].get_sublayer_boundaries(prev_depth, bottom)
        prev_depth = bottom
    return boundaries

get_sublayer_depths()

Get a sorted list of all sublayer boundary depths in the profile.

Source code in simplesoilprofile/models/profile.py 
80
81
82
83
84
85
def get_sublayer_depths(self) -> list[float]:
    """Get a sorted list of all sublayer boundary depths in the profile."""
    all_depths = []
    for depths in self.get_sublayer_boundaries().values():
        all_depths.extend(depths)
    return sorted(set(all_depths))

validate_layer_depths()

Validate that layer depths are consistent and monotonically increasing.

Source code in simplesoilprofile/models/profile.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
@model_validator(mode='after')
def validate_layer_depths(self) -> 'SoilProfile':
    """Validate that layer depths are consistent and monotonically increasing."""
    if not self.layer_bottoms:
        raise ValueError("Layer depths cannot be empty")

    if len(self.layer_bottoms) != len(self.layers):
        raise ValueError("Must provide bottom depths for all layers")

    prev_depth = 0
    for i, bottom in enumerate(self.layer_bottoms):
        if bottom <= prev_depth:
            raise ValueError(f"Layer {i}: bottom depth must be greater than previous layer")
        prev_depth = bottom

    return self

get_profile_from_dov(location, fetch_elevation=False, crs='EPSG:31370')

Fetch a soil profile from DOV WMS at a specific location.

This function queries the DOV WMS service for soil texture data at the given location and constructs a SoilProfile object.

Parameters:

Name Type Description Default
location Point

Shapely Point representing the location (x, y coordinates)

 required
fetch_elevation bool

Whether to fetch elevation data for the profile surface

 False
crs str

Coordinate reference system of the input location

 'EPSG:31370'

Returns:

Type Description
SoilProfile | None

SoilProfile object or None if data not found
Source code in simplesoilprofile/models/profile.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def get_profile_from_dov(
    location: Point,
    fetch_elevation: bool = False,
    crs: str = "EPSG:31370"
) -> SoilProfile | None:
    """Fetch a soil profile from DOV WMS at a specific location.

    This function queries the DOV WMS service for soil texture data at the given
    location and constructs a SoilProfile object.

    Args:
        location: Shapely Point representing the location (x, y coordinates)
        fetch_elevation: Whether to fetch elevation data for the profile surface
        crs: Coordinate reference system of the input location

    Returns:
        SoilProfile object or None if data not found
    """
    try:
        from dovwms import get_profile_from_dov
    except ImportError as e:
        raise ImportError("Please install dovwms using 'pip install dovwms'.") from e

    profile_data = get_profile_from_dov(
        location.x,
        location.y,
        fetch_elevation=fetch_elevation,
        crs=crs
    )

    profile = SoilProfile(
        name="DOV Soil Profile",
        location=location,
        elevation=profile_data.get("elevation").get("elevation") if fetch_elevation else None,
        layers=[
            SoilLayer(
                name=layer["name"],
                layer_top=layer["layer_top"],
                layer_bottom=layer["layer_bottom"],
                sand_content=layer["sand_content"],
                silt_content=layer["silt_content"],
                clay_content=layer["clay_content"],
                metadata=layer.get("metadata", {})
            )
            for layer in profile_data["layers"]
        ],
        layer_bottoms=[layer["layer_bottom"] for layer in profile_data["layers"]]
    )
    return profile

Soil layer model representing a single layer with physical properties.

SoilLayer

Bases: BaseModel

A soil layer with uniform properties.

This class represents a soil layer with van Genuchten parameters and other properties needed for hydrological modeling, particularly for SWAP model integration. The layer can be discretized into sublayers for numerical computation.

Source code in simplesoilprofile/models/layer.py 
 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
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
class SoilLayer(BaseModel):
    """A soil layer with uniform properties.

    This class represents a soil layer with van Genuchten parameters and other
    properties needed for hydrological modeling, particularly for SWAP model
    integration. The layer can be discretized into sublayers for numerical computation.
    """

    name: str = Field(..., description="Name or identifier of the soil layer")
    texture_class: str | None = Field(None, description="Soil texture class (e.g., 'sand', 'clay', etc.)")
    discretization: LayerDiscretization | None = Field(
        None, description="Configuration for layer discretization into sublayers"
    )
    description: str | None = Field(None, description="Description of the soil layer")

    # Van Genuchten parameters
    theta_res: float | None = Field(None, description="Residual water content [cm³/cm³]", ge=0.0, le=1.0)
    theta_sat: float | None = Field(None, description="Saturated water content [cm³/cm³]", ge=0.0, le=1.0)
    alpha: float | None = Field(None, description="Shape parameter alpha [1/cm]", gt=0.0)
    n: float | None = Field(None, description="Shape parameter n [-]", gt=1.0)
    k_sat: float | None = Field(None, description="Saturated hydraulic conductivity [cm/day]", gt=0.0)
    lambda_param: float | None = Field(0.5, description="Tortuosity parameter lambda [-]", alias="l")
    alphaw: float | None = Field(None, description="Alfa parameter of main wetting curve [cm]", ge=0.0, le=100.0)
    h_enpr: float | None = Field(None, description="Air entry pressure head [cm]", ge=-40.0, le=0)
    ksatexm: float | None = Field(None, description="Measured hydraulic conductivity at saturated conditions [cm/day]", ge=0.0, le=1e5)
    bulk_density: float | None = Field(None, description="Bulk density [g/cm³]", gt=0.0)

    #  texture information (used in heat flow section of SWAP)
    clay_content: float | None = Field(None, description="Clay content [%]", ge=0.0, le=100.0)
    silt_content: float | None = Field(None, description="Silt content [%]", ge=0.0, le=100.0)
    sand_content: float | None = Field(None, description="Sand content [%]", ge=0.0, le=100.0)
    organic_matter: float | None = Field(None, description="Organic matter content [%]", ge=0.0)

    metadata: dict[str, M] | None = Field(
        default_factory=dict,
        description="Metadata for each measured parameter, keyed by parameter name"
    )

    @model_validator(mode='after')
    def validate_water_contents(self) -> 'SoilLayer':
        """Validate that residual water content is less than saturated water content."""
        if not self.theta_res or not self.theta_sat:
            return self

        if self.theta_res >= self.theta_sat:
            raise ValueError("Residual water content must be less than saturated water content")
        return self

    @computed_field
    @property
    def sum_texture(self) -> float | None:
        """Compute the sum of clay, silt, and sand contents if all are provided."""
        if all(v is not None for v in (self.clay_content, self.silt_content, self.sand_content)):
            return self.clay_content + self.silt_content + self.sand_content
        return None

    def normalize_soil_fractions(
        self,
        tolerance: float = 2.0
    ) -> tuple[float, float, float]:
        """
        Normalize soil texture fractions to sum to 100%.

        Parameters
        ----------
        sand, silt, clay : float
            Soil fraction percentages
        tolerance : float
            Maximum acceptable deviation from 100% before normalization (default: 2%)

        Returns
        -------
        tuple[float, float, float]
            Normalized (sand, silt, clay) percentages

        Examples
        --------
        >>> normalize_soil_fractions(72.97, 22.44, 3.61)
        (73.76, 22.68, 3.65)  # Sum was 99.02, now 100.0
        """
        total = self.sand_content + self.silt_content + self.clay_content

        # Check if normalization is needed
        if abs(total - 100.0) < 0.01:  # Already at 100%
            return None

        # Warn if deviation is large (potential data quality issue)
        if abs(total - 100.0) > tolerance:
            logger.warning(
                "Sum of soil fractions %.2f%% exceeds tolerance of %.1f%% for layer '%s'",
                total, tolerance, self.name
            )

        # Proportional normalization
        if total > 0:
            logger.warning(
                "Normalizing soil fractions for layer '%s': %.2f%% -> 100.0%%",
                self.name, total
            )
            self.sand_content = (self.sand_content / total) * 100.0
            self.silt_content = (self.silt_content / total) * 100.0
            self.clay_content = (self.clay_content / total) * 100.0

            return None
        else:
            raise ValueError("Sum of fractions is zero or negative")

    def get_sublayer_boundaries(self, top: float, bottom: float) -> list[float]:
        """Get the sublayer boundary depths for this layer.

        Args:
            top: Top depth of the layer [cm]
            bottom: Bottom depth of the layer [cm]

        Returns:
            List of boundary depths [cm], including top and bottom depths.
            If no discretization is configured, returns [top, bottom].
        """
        if self.discretization is None:
            return [top, bottom]

        return compute_sublayer_boundaries(top, bottom, self.discretization)

    def predict_van_genuchten(self, method: Literal["rosetta",]):
        """Predict van Genuchten parameters from soil texture."""
        if method == "rosetta":
            input_data = [
                [self.sand_content, self.silt_content, self.clay_content]
            ]
            soildata = SoilData.from_array(input_data)

            mean, _stdev, _codes = rosetta(2, soildata)
            logger.debug("Raw Rosetta output (mean values): %s", mean)

            # Convert from log10 values for alpha, n, and k_sat
            self.theta_res = mean[0][0]  # Residual water content
            self.theta_sat = mean[0][1]  # Saturated water content
            self.alpha = 10 ** mean[0][2]  # Convert from log10(1/cm)
            self.n = 10 ** mean[0][3]      # Convert from log10(-)
            self.k_sat = 10 ** mean[0][4]  # Convert from log10(cm/day)

            logger.info("Predicted van Genuchten parameters for layer '%s'")

    def infer_fractions_from_texture(
            self,
            texture_class: str
        ) -> None:
        """
        Add sand/silt/clay percentages to a SoilLayer based on texture class.

        Parameters
        ----------
        layer : SoilLayer
            Soil layer to enrich
        texture_class : str
            Texture class name from field observations

        Returns
        -------
        SoilLayer
            Updated layer with estimated percentages
        """
        project_root = os.path.dirname(os.path.dirname(os.path.dirname(__file__)))
        data_path = os.path.join(project_root, 'simplesoilprofile', 'models', 'data', 'usda_texture.yaml')
        texture_converter = SoilTextureConverter(data_path)

        # Get percentages from texture class
        sand, silt, clay = texture_converter.class_to_percentages(texture_class)

        # Update layer
        self.sand_content = sand
        self.silt_content = silt
        self.clay_content = clay
        self.texture_class = texture_class

        # Add metadata about the conversion
        self.description = (
            f"{self.description or ''} "
            f"[Texture fractions estimated from class '{texture_class}' "
            f"using {texture_converter.metadata['reference']}]"
        )

        return None

sum_texture property

Compute the sum of clay, silt, and sand contents if all are provided.

get_sublayer_boundaries(top, bottom)

Get the sublayer boundary depths for this layer.

Parameters:

Name Type Description Default
top float

Top depth of the layer [cm]

 required
bottom float

Bottom depth of the layer [cm]

 required

Returns:

Type Description
list[float]

List of boundary depths [cm], including top and bottom depths.
list[float]

If no discretization is configured, returns [top, bottom].
Source code in simplesoilprofile/models/layer.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def get_sublayer_boundaries(self, top: float, bottom: float) -> list[float]:
    """Get the sublayer boundary depths for this layer.

    Args:
        top: Top depth of the layer [cm]
        bottom: Bottom depth of the layer [cm]

    Returns:
        List of boundary depths [cm], including top and bottom depths.
        If no discretization is configured, returns [top, bottom].
    """
    if self.discretization is None:
        return [top, bottom]

    return compute_sublayer_boundaries(top, bottom, self.discretization)

infer_fractions_from_texture(texture_class)

Add sand/silt/clay percentages to a SoilLayer based on texture class.

Parameters

layer : SoilLayer Soil layer to enrich texture_class : str Texture class name from field observations

Returns

SoilLayer Updated layer with estimated percentages

Source code in simplesoilprofile/models/layer.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
def infer_fractions_from_texture(
        self,
        texture_class: str
    ) -> None:
    """
    Add sand/silt/clay percentages to a SoilLayer based on texture class.

    Parameters
    ----------
    layer : SoilLayer
        Soil layer to enrich
    texture_class : str
        Texture class name from field observations

    Returns
    -------
    SoilLayer
        Updated layer with estimated percentages
    """
    project_root = os.path.dirname(os.path.dirname(os.path.dirname(__file__)))
    data_path = os.path.join(project_root, 'simplesoilprofile', 'models', 'data', 'usda_texture.yaml')
    texture_converter = SoilTextureConverter(data_path)

    # Get percentages from texture class
    sand, silt, clay = texture_converter.class_to_percentages(texture_class)

    # Update layer
    self.sand_content = sand
    self.silt_content = silt
    self.clay_content = clay
    self.texture_class = texture_class

    # Add metadata about the conversion
    self.description = (
        f"{self.description or ''} "
        f"[Texture fractions estimated from class '{texture_class}' "
        f"using {texture_converter.metadata['reference']}]"
    )

    return None

normalize_soil_fractions(tolerance=2.0)

Normalize soil texture fractions to sum to 100%.

Parameters

sand, silt, clay : float Soil fraction percentages tolerance : float Maximum acceptable deviation from 100% before normalization (default: 2%)

Returns

tuple[float, float, float] Normalized (sand, silt, clay) percentages

Examples

normalize_soil_fractions(72.97, 22.44, 3.61) (73.76, 22.68, 3.65) # Sum was 99.02, now 100.0

Source code in simplesoilprofile/models/layer.py 
 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
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def normalize_soil_fractions(
    self,
    tolerance: float = 2.0
) -> tuple[float, float, float]:
    """
    Normalize soil texture fractions to sum to 100%.

    Parameters
    ----------
    sand, silt, clay : float
        Soil fraction percentages
    tolerance : float
        Maximum acceptable deviation from 100% before normalization (default: 2%)

    Returns
    -------
    tuple[float, float, float]
        Normalized (sand, silt, clay) percentages

    Examples
    --------
    >>> normalize_soil_fractions(72.97, 22.44, 3.61)
    (73.76, 22.68, 3.65)  # Sum was 99.02, now 100.0
    """
    total = self.sand_content + self.silt_content + self.clay_content

    # Check if normalization is needed
    if abs(total - 100.0) < 0.01:  # Already at 100%
        return None

    # Warn if deviation is large (potential data quality issue)
    if abs(total - 100.0) > tolerance:
        logger.warning(
            "Sum of soil fractions %.2f%% exceeds tolerance of %.1f%% for layer '%s'",
            total, tolerance, self.name
        )

    # Proportional normalization
    if total > 0:
        logger.warning(
            "Normalizing soil fractions for layer '%s': %.2f%% -> 100.0%%",
            self.name, total
        )
        self.sand_content = (self.sand_content / total) * 100.0
        self.silt_content = (self.silt_content / total) * 100.0
        self.clay_content = (self.clay_content / total) * 100.0

        return None
    else:
        raise ValueError("Sum of fractions is zero or negative")

predict_van_genuchten(method)

Predict van Genuchten parameters from soil texture.

Source code in simplesoilprofile/models/layer.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def predict_van_genuchten(self, method: Literal["rosetta",]):
    """Predict van Genuchten parameters from soil texture."""
    if method == "rosetta":
        input_data = [
            [self.sand_content, self.silt_content, self.clay_content]
        ]
        soildata = SoilData.from_array(input_data)

        mean, _stdev, _codes = rosetta(2, soildata)
        logger.debug("Raw Rosetta output (mean values): %s", mean)

        # Convert from log10 values for alpha, n, and k_sat
        self.theta_res = mean[0][0]  # Residual water content
        self.theta_sat = mean[0][1]  # Saturated water content
        self.alpha = 10 ** mean[0][2]  # Convert from log10(1/cm)
        self.n = 10 ** mean[0][3]      # Convert from log10(-)
        self.k_sat = 10 ** mean[0][4]  # Convert from log10(cm/day)

        logger.info("Predicted van Genuchten parameters for layer '%s'")

validate_water_contents()

Validate that residual water content is less than saturated water content.

Source code in simplesoilprofile/models/layer.py 
55
56
57
58
59
60
61
62
63
@model_validator(mode='after')
def validate_water_contents(self) -> 'SoilLayer':
    """Validate that residual water content is less than saturated water content."""
    if not self.theta_res or not self.theta_sat:
        return self

    if self.theta_res >= self.theta_sat:
        raise ValueError("Residual water content must be less than saturated water content")
    return self

Layer discretization models and utilities.

DiscretizationType

Bases: str, Enum

Types of layer discretization methods.

Source code in simplesoilprofile/models/discretization.py 
 9
10
11
12
13
14
15
class DiscretizationType(str, Enum):
    """Types of layer discretization methods."""

    EVEN = "even"
    LOG_TOP = "logarithmic_top"
    LOG_BOTTOM = "logarithmic_bottom"
    LOG_BOTH = "logarithmic_both"

LayerDiscretization

Bases: BaseModel

Configuration for layer discretization.

Source code in simplesoilprofile/models/discretization.py 
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
class LayerDiscretization(BaseModel):
    """Configuration for layer discretization."""

    type: DiscretizationType = Field(
        ...,
        description="Type of discretization to apply"
    )
    num_sublayers: int = Field(
        ...,
        description="Number of sublayers to create",
        gt=0
    )
    num_compartments: int = Field(
        ...,
        description="Number of compartments to create",
        gt=0
    )
    log_density: float = Field(
        1.0,
        description="Density parameter for logarithmic spacing (higher = more dense)",
        gt=0
    )

    @model_validator(mode='after')
    def validate_discretization(self) -> 'LayerDiscretization':
        """Validate discretization parameters."""
        if self.type != DiscretizationType.EVEN and self.log_density <= 0:
            raise ValueError("log_density must be positive for logarithmic discretization")
        return self

    @computed_field
    @property
    def compartment_heights(self) -> list[float]:
        """Get the ordered list of compartment heights for this discretization.

        The heights are calculated based on the discretization type and parameters.
        The order is from top to bottom of the layer.

        Returns:
            List of compartment heights [cm] in order from top to bottom.
            For EVEN discretization, all heights are equal.
            For logarithmic discretizations, heights vary according to the spacing function.
        """
        # Get sublayer boundaries for a unit thickness (1cm) layer
        # (we'll use the ratios, actual heights will be scaled by layer thickness)
        boundaries = compute_sublayer_boundaries(0, 1, self)
        heights = list(np.diff(boundaries))
        return heights

compartment_heights property

Get the ordered list of compartment heights for this discretization.

The heights are calculated based on the discretization type and parameters. The order is from top to bottom of the layer.

Returns:

Type Description
list[float]

List of compartment heights [cm] in order from top to bottom.
list[float]

For EVEN discretization, all heights are equal.
list[float]

For logarithmic discretizations, heights vary according to the spacing function.

validate_discretization()

Validate discretization parameters.

Source code in simplesoilprofile/models/discretization.py 
40
41
42
43
44
45
@model_validator(mode='after')
def validate_discretization(self) -> 'LayerDiscretization':
    """Validate discretization parameters."""
    if self.type != DiscretizationType.EVEN and self.log_density <= 0:
        raise ValueError("log_density must be positive for logarithmic discretization")
    return self

compute_sublayer_boundaries(top, bottom, discretization)

Compute sublayer boundaries based on discretization configuration.

Parameters:

Name Type Description Default
top float

Top depth of the layer [cm]

 required
bottom float

Bottom depth of the layer [cm]

 required
discretization LayerDiscretization

Discretization configuration

 required

Returns:

Type Description
list[float]

List of boundary depths [cm], including top and bottom
Source code in simplesoilprofile/models/discretization.py 
 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
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def compute_sublayer_boundaries(
    top: float,
    bottom: float,
    discretization: LayerDiscretization
) -> list[float]:
    """Compute sublayer boundaries based on discretization configuration.

    Args:
        top: Top depth of the layer [cm]
        bottom: Bottom depth of the layer [cm]
        discretization: Discretization configuration

    Returns:
        List of boundary depths [cm], including top and bottom
    """
    thickness = bottom - top
    n = discretization.num_sublayers

    if discretization.type == DiscretizationType.EVEN:
        # Simple linear spacing
        return list(np.linspace(top, bottom, n + 1))

    elif discretization.type == DiscretizationType.LOG_TOP:
        # Logarithmic spacing with finer discretization at the top
        # Use exponential function to create non-uniform spacing
        positions = np.exp(np.linspace(0, discretization.log_density, n + 1)) - 1
        positions = positions / positions[-1]  # normalize to [0, 1]
        return list(top + positions * thickness)

    elif discretization.type == DiscretizationType.LOG_BOTTOM:
        # Logarithmic spacing with finer discretization at the bottom
        # Use reversed exponential function
        positions = np.exp(np.linspace(0, discretization.log_density, n + 1)) - 1
        positions = positions / positions[-1]  # normalize to [0, 1]
        positions = 1 - positions[::-1]  # reverse and invert
        return list(top + positions * thickness)

    else:  # LOG_BOTH
        # Symmetric logarithmic spacing with finer discretization at both ends
        if n % 2 == 0:
            n += 1  # ensure odd number for middle point

        # Split the layer into two parts
        middle = (top + bottom) / 2
        n_half = (n + 1) // 2

        # Create top half (fine → coarse)
        top_positions = np.exp(np.linspace(0, discretization.log_density, n_half)) - 1
        top_positions = top_positions / top_positions[-1]
        top_boundaries = top + top_positions * (middle - top)

        # Create bottom half (coarse → fine)
        bottom_positions = 1 - top_positions[::-1]
        bottom_boundaries = middle + bottom_positions * (bottom - middle)

        # Combine boundaries, avoiding duplicate middle point
        return list(np.concatenate([top_boundaries[:-1], bottom_boundaries]))

SoilLayerMetadata

Bases: BaseModel

Metadata describing the source and quality of a measurement.

This class captures provenance information for soil measurements, following best practices for scientific data traceability and reproducibility.

Source code in simplesoilprofile/models/metadata.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
class SoilLayerMetadata(BaseModel):
    """Metadata describing the source and quality of a measurement.

    This class captures provenance information for soil measurements,
    following best practices for scientific data traceability and reproducibility.
    """

    # Data source identification
    source: str | None = Field("User-provided", description="Source of the data (e.g., 'DOV API', 'Field measurement', 'Laboratory analysis')")
    url: HttpUrl | None = Field(None, description="URL to original data source")
    source_type: Literal['measured', 'modeled', 'derived', 'literature', 'estimated', "unknown"] | None = Field(
        "unknown", description="Type of data source"
    )
    access_date: date | None = Field(None, description="Date when data was retrieved/accessed")
    uncertainty: float | None = Field(None, description="Measurement uncertainty or confidence interval")

Plotting

Plotting utilities for visualizing soil profiles.

plot_profile(profile, ax=None, figsize=(8, 12), texture_colors=None, show_depths=True, show_layer_properties=True, show_sublayers=True)

Plot a soil profile showing layers and their properties.

Parameters:

Name Type Description Default
profile SoilProfile

The soil profile to plot

 required
ax Axes | None

Optional matplotlib axes to plot on

 None
figsize tuple[float, float]

Figure size (width, height) if creating new figure

 (8, 12)
texture_colors dict[str, str] | None

Optional mapping of texture classes to colors

 None
show_depths bool

Whether to show depth labels

 True
show_properties

Whether to show soil properties in annotations

 required
show_sublayers bool

Whether to show sublayer boundaries for discretized layers

 True

Returns:

Type Description
Axes

The matplotlib axes object containing the plot
Source code in simplesoilprofile/plotting/profile_plot.py 
 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
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def plot_profile(
    profile: SoilProfile,
    ax: plt.Axes | None = None,
    figsize: tuple[float, float] = (8, 12),
    texture_colors: dict[str, str] | None = None,
    show_depths: bool = True,
    show_layer_properties: bool = True,
    show_sublayers: bool = True,
) -> plt.Axes:
    """Plot a soil profile showing layers and their properties.

    Args:
        profile: The soil profile to plot
        ax: Optional matplotlib axes to plot on
        figsize: Figure size (width, height) if creating new figure
        texture_colors: Optional mapping of texture classes to colors
        show_depths: Whether to show depth labels
        show_properties: Whether to show soil properties in annotations
        show_sublayers: Whether to show sublayer boundaries for discretized layers

    Returns:
        The matplotlib axes object containing the plot
    """

    def _show_depths(ax: plt.Axes):
        ax.yaxis.set_major_locator(plt.MultipleLocator(20))
        ax.yaxis.set_minor_locator(plt.MultipleLocator(5))
        ax.grid(True, axis='y', which='major', linestyle='--', alpha=0.3)

    def _format_layer_properties(layer: "SoilLayer", ax: plt.Axes):
        props = [
            f"Layer: {layer.name}",
            f"Texture: {layer.texture_class or 'Unknown'}"
        ]

        if all(x is not None for x in [layer.clay_content, layer.silt_content, layer.sand_content]):
            props.append(f"Clay/Silt/Sand: {layer.clay_content:.0f}/{layer.silt_content:.0f}/{layer.sand_content:.0f}%")

        # Only show hydraulic properties if they exist
        if layer.theta_res is not None and layer.theta_sat is not None:
            props.append(f"θr/θs: {layer.theta_res:.3f}/{layer.theta_sat:.3f}")

        if layer.k_sat is not None:
            props.append(f"Ks: {layer.k_sat:.2f} cm/d")

        # Add text annotation
        mid_depth = (top + bottom) / 2
        ax.text(
            1.1, -mid_depth,
            '\n'.join(props),
            verticalalignment='center',
            horizontalalignment='left',
            fontsize=8,
        )

    def _show_sublayers(layer: "SoilLayer", top: float, bottom: float, ax: plt.Axes):
        sublayer_depths = layer.get_sublayer_boundaries(top, bottom)
        for depth in sublayer_depths[1:-1]:  # Skip top and bottom depths
            ax.axhline(
                y=-depth,
                xmin=0,
                xmax=1,
                color='orange',
                linestyle=':',
                linewidth=0.5,
                alpha=0.5
            )

    if ax is None:
        _, ax = plt.subplots(figsize=figsize)

    if texture_colors is None:
        texture_colors = DEFAULT_TEXTURE_COLORS

    # Sort layers by depth
    sorted_depths = sorted(profile.layer_bounds.items(), key=lambda x: x[0])
    total_depth = profile.profile_depth

    # Plot each layer
    for layer_idx, (top, bottom) in sorted_depths:
        layer = profile.layers[layer_idx]

        # Get color based on texture class or use default
        color = texture_colors.get(
            layer.texture_class.lower() if layer.texture_class else 'unknown',
            '#808080'  # Default gray for unknown textures
        )

        # Plot the layer as a rectangle
        rect = plt.Rectangle(
            (0, -bottom),           # (x, y) of bottom-left corner
            1,                      # width (normalized to 1)
            bottom - top,           # height
            facecolor=color,
            edgecolor='black',
            linewidth=0.5,
        )
        ax.add_patch(rect)

        # Add sublayer lines if enabled and layer has discretization
        if show_sublayers and layer.discretization is not None:
            _show_sublayers(layer, top, bottom, ax)

        # Add layer information if requested
        if show_layer_properties:
            _format_layer_properties(layer, ax)

    # Set axis limits and labels
    ax.set_xlim(-0.1, 2.0)  # Leave space for annotations
    ax.set_ylim(-total_depth * 1.1, 0)  # Add 10% padding at bottom

    # Remove x-axis and right spine
    ax.xaxis.set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_visible(False)

    # Add depth ticks if requested
    if show_depths:
        _show_depths(ax=ax)


    # Add title with profile information
    title = f"Soil Profile: {profile.name}"
    if profile.location is not None and any(coord is not None for coord in [profile.location.x, profile.location.y, profile.elevation]):
        coords = []
        if profile.location.x is not None and profile.location.y is not None:
            coords.append(f"x={profile.location.x:.1f}, y={profile.location.y:.1f}")
        if profile.elevation is not None:
            coords.append(f"z={profile.elevation:.1f}")
        if coords:
            title += f"\n({', '.join(coords)})"

    ax.set_title(title)
    ax.set_ylabel('Depth [cm]')

    return ax

SWAP Integration

Utilities for generating SWAP model input files from soil profiles.

profile_to_soilhydfunc_table(profile)

Convert a SoilProfile to a SWAP-compatible SOILHYDRFUNC table.

Parameters:

Name Type Description Default
profile SoilProfile

The soil profile to convert

 required

Returns:

Type Description
DataFrame

DataFrame containing the SOILHYDRFUNC table with any column that contains
DataFrame

None/NaN values removed.
Source code in simplesoilprofile/models/swap.py 
 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
def profile_to_soilhydfunc_table(
    profile: SoilProfile,
) -> pd.DataFrame:
    """Convert a SoilProfile to a SWAP-compatible SOILHYDRFUNC table.

    Args:
        profile: The soil profile to convert

    Returns:
        DataFrame containing the SOILHYDRFUNC table with any column that contains
        None/NaN values removed.
    """
    rows = []
    for layer in profile.layers:
        rows.append({
            "ORES": layer.theta_res,
            "OSAT": layer.theta_sat,
            "ALFA": layer.alpha,
            "NPAR": layer.n,
            "LEXP": layer.lambda_param,
            "KSATFIT": layer.k_sat,
            "H_ENPR": layer.h_enpr,
            "KSATEXM": layer.ksatexm,
            "BDENS": layer.bulk_density,
            "ALFAW": layer.alphaw,
        })

    df = pd.DataFrame(rows)
    df = df.dropna(axis=1, how="any")
    return df

profile_to_sublayer_table(profile)

Convert a SoilProfile to a DataFrame representing sublayers and compartments.

Source code in simplesoilprofile/models/swap.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def profile_to_sublayer_table(profile: SoilProfile) -> pd.DataFrame:
    """Convert a SoilProfile to a DataFrame representing sublayers and compartments."""
    rows = []
    sublay_counter = 1
    for isoillay, layer in enumerate(profile.layers, start=1):
        if layer.discretization:
            # scale normalized compartment heights by the actual layer thickness
            top, bottom = profile.layer_bounds[isoillay - 1]
            thickness = bottom - top
            for height_ratio in layer.discretization.compartment_heights:
                height = height_ratio * thickness
                rows.append({
                    'ISUBLAY': sublay_counter,
                    'ISOILLAY': isoillay,
                    'HSUBLAY': height,
                    'NCOMP': layer.discretization.num_compartments,
                    'HCOMP': height / layer.discretization.num_compartments
                })
                sublay_counter += 1

    return pd.DataFrame(rows)

profile_to_texture_table(profile)

Convert a SoilProfile to a SWAP-compatible SOILTEXTURE table.

Parameters:

Name Type Description Default
profile SoilProfile

The soil profile to convert

 required

Returns:

Type Description
DataFrame

DataFrame containing the SOILTEXTURE table.
Source code in simplesoilprofile/models/swap.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def profile_to_texture_table(
    profile: SoilProfile,
) -> pd.DataFrame:
    """Convert a SoilProfile to a SWAP-compatible SOILTEXTURE table.

    Args:
        profile: The soil profile to convert

    Returns:
        DataFrame containing the SOILTEXTURE table.
    """
    rows = []
    for layer in profile.layers:
        rows.append({
            "PCLAY": layer.clay_content,
            "PSILT": layer.silt_content,
            "PSAND": layer.sand_content,
            "ORGMAT": layer.organic_matter,
        })
    return pd.DataFrame(rows)

Texture Conversion

SoilTextureConverter

Convert between soil texture class names and sand/silt/clay percentages.

Uses USDA classification data loaded from YAML configuration file.

Source code in simplesoilprofile/models/texture_conversion.py 
  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
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
class SoilTextureConverter:
    """Convert between soil texture class names and sand/silt/clay percentages.

    Uses USDA classification data loaded from YAML configuration file.
    """

    def __init__(self, config_path: str = 'usda_texture_data.yaml'):
        """
        Initialize with YAML configuration file.

        Parameters
        ----------
        config_path : str
            Path to YAML configuration file
        """
        with open(config_path) as f:
            self.config = yaml.safe_load(f)

        self.centroids = self.config['centroids']
        self.ranges = self.config['ranges']
        self.aliases = self.config['aliases']
        self.metadata = self.config['metadata']

    def _normalize_class_name(self, texture_class: str) -> str:
        """Normalize texture class name to match YAML keys."""
        # Convert to lowercase and replace spaces with underscores
        normalized = texture_class.lower().strip().replace(' ', '_').replace('-', '_')

        # Check if it's an alias
        if normalized in self.aliases:
            return self.aliases[normalized]

        return normalized

    def class_to_percentages(
        self,
        texture_class: str,
        method: str = 'centroid',
        normalize: bool = True
    ) -> tuple[float, float, float]:
        """
        Convert texture class name to sand/silt/clay percentages.

        Parameters
        ----------
        texture_class : str
            Soil texture class name (e.g., 'silty sand', 'sandy loam')
        method : str
            'centroid' for geometric centroid or 'mean' for statistical mean
        normalize : bool
            Whether to normalize to sum to 100%

        Returns
        -------
        tuple[float, float, float]
            (sand, silt, clay) percentages

        Examples
        --------
        >>> converter = SoilTextureConverter()
        >>> sand, silt, clay = converter.class_to_percentages('loamy sand')
        >>> print(f"Sand: {sand}%, Silt: {silt}%, Clay: {clay}%")
        Sand: 82.0%, Silt: 12.0%, Clay: 6.0%
        """
        # Normalize the class name
        class_key = self._normalize_class_name(texture_class)

        # Get data based on method
        if method == 'centroid':
            if class_key not in self.centroids:
                raise ValueError(
                    f"Unknown texture class: '{texture_class}'. "
                    f"Available classes: {list(self.centroids.keys())}"
                )
            data = self.centroids[class_key]
            sand = data['sand']
            silt = data['silt']
            clay = data['clay']

        elif method == 'mean':
            if class_key not in self.ranges:
                raise ValueError(
                    f"Unknown texture class: '{texture_class}'. "
                    f"Available classes: {list(self.ranges.keys())}"
                )
            data = self.ranges[class_key]
            sand = data['sand']['mean']
            silt = data['silt']['mean']
            clay = data['clay']['mean']

        else:
            raise ValueError(f"Unknown method: '{method}'. Use 'centroid' or 'mean'")

        # Normalize if requested
        if normalize:
            total = sand + silt + clay
            if total > 0:
                sand = (sand / total) * 100
                silt = (silt / total) * 100
                clay = (clay / total) * 100

        return sand, silt, clay

    def get_ranges(self, texture_class: str) -> dict:
        """
        Get statistical ranges for a texture class.

        Parameters
        ----------
        texture_class : str
            Soil texture class name

        Returns
        -------
        dict
            Dictionary with mean, min, max, std for each fraction
        """
        class_key = self._normalize_class_name(texture_class)

        if class_key not in self.ranges:
            raise ValueError(f"Unknown texture class: '{texture_class}'")

        return self.ranges[class_key]

    def get_metadata(self) -> dict:
        """Get metadata about the classification system."""
        return self.metadata

    def list_available_classes(self) -> list:
        """List all available texture classes."""
        return list(self.centroids.keys())

__init__(config_path='usda_texture_data.yaml')

Initialize with YAML configuration file.

Parameters

config_path : str Path to YAML configuration file

Source code in simplesoilprofile/models/texture_conversion.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def __init__(self, config_path: str = 'usda_texture_data.yaml'):
    """
    Initialize with YAML configuration file.

    Parameters
    ----------
    config_path : str
        Path to YAML configuration file
    """
    with open(config_path) as f:
        self.config = yaml.safe_load(f)

    self.centroids = self.config['centroids']
    self.ranges = self.config['ranges']
    self.aliases = self.config['aliases']
    self.metadata = self.config['metadata']

class_to_percentages(texture_class, method='centroid', normalize=True)

Convert texture class name to sand/silt/clay percentages.

Parameters

texture_class : str Soil texture class name (e.g., 'silty sand', 'sandy loam') method : str 'centroid' for geometric centroid or 'mean' for statistical mean normalize : bool Whether to normalize to sum to 100%

Returns

tuple[float, float, float] (sand, silt, clay) percentages

Examples

converter = SoilTextureConverter() sand, silt, clay = converter.class_to_percentages('loamy sand') print(f"Sand: {sand}%, Silt: {silt}%, Clay: {clay}%") Sand: 82.0%, Silt: 12.0%, Clay: 6.0%

Source code in simplesoilprofile/models/texture_conversion.py 
 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
 99
100
101
102
103
104
105
106
def class_to_percentages(
    self,
    texture_class: str,
    method: str = 'centroid',
    normalize: bool = True
) -> tuple[float, float, float]:
    """
    Convert texture class name to sand/silt/clay percentages.

    Parameters
    ----------
    texture_class : str
        Soil texture class name (e.g., 'silty sand', 'sandy loam')
    method : str
        'centroid' for geometric centroid or 'mean' for statistical mean
    normalize : bool
        Whether to normalize to sum to 100%

    Returns
    -------
    tuple[float, float, float]
        (sand, silt, clay) percentages

    Examples
    --------
    >>> converter = SoilTextureConverter()
    >>> sand, silt, clay = converter.class_to_percentages('loamy sand')
    >>> print(f"Sand: {sand}%, Silt: {silt}%, Clay: {clay}%")
    Sand: 82.0%, Silt: 12.0%, Clay: 6.0%
    """
    # Normalize the class name
    class_key = self._normalize_class_name(texture_class)

    # Get data based on method
    if method == 'centroid':
        if class_key not in self.centroids:
            raise ValueError(
                f"Unknown texture class: '{texture_class}'. "
                f"Available classes: {list(self.centroids.keys())}"
            )
        data = self.centroids[class_key]
        sand = data['sand']
        silt = data['silt']
        clay = data['clay']

    elif method == 'mean':
        if class_key not in self.ranges:
            raise ValueError(
                f"Unknown texture class: '{texture_class}'. "
                f"Available classes: {list(self.ranges.keys())}"
            )
        data = self.ranges[class_key]
        sand = data['sand']['mean']
        silt = data['silt']['mean']
        clay = data['clay']['mean']

    else:
        raise ValueError(f"Unknown method: '{method}'. Use 'centroid' or 'mean'")

    # Normalize if requested
    if normalize:
        total = sand + silt + clay
        if total > 0:
            sand = (sand / total) * 100
            silt = (silt / total) * 100
            clay = (clay / total) * 100

    return sand, silt, clay

get_metadata()

Get metadata about the classification system.

Source code in simplesoilprofile/models/texture_conversion.py 
129
130
131
def get_metadata(self) -> dict:
    """Get metadata about the classification system."""
    return self.metadata

get_ranges(texture_class)

Get statistical ranges for a texture class.

Parameters

texture_class : str Soil texture class name

Returns

dict Dictionary with mean, min, max, std for each fraction

Source code in simplesoilprofile/models/texture_conversion.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def get_ranges(self, texture_class: str) -> dict:
    """
    Get statistical ranges for a texture class.

    Parameters
    ----------
    texture_class : str
        Soil texture class name

    Returns
    -------
    dict
        Dictionary with mean, min, max, std for each fraction
    """
    class_key = self._normalize_class_name(texture_class)

    if class_key not in self.ranges:
        raise ValueError(f"Unknown texture class: '{texture_class}'")

    return self.ranges[class_key]

list_available_classes()

List all available texture classes.

Source code in simplesoilprofile/models/texture_conversion.py 
133
134
135
def list_available_classes(self) -> list:
    """List all available texture classes."""
    return list(self.centroids.keys())

Utilities

Utility functions and configuration for the simplesoilprofile package.