Soil water-holding capacity - gridded for MÉRA

European soil database derived data - total available water content (TAWC) for the topsoil [mm] (European Commission, n.d.; Hiederer, 2013a; Hiederer, 2013b): https://esdac.jrc.ec.europa.eu/content/european-soil-database-derived-data

import os
import matplotlib.pyplot as plt
import geopandas as gpd
import rioxarray as rxr
from rasterstats import zonal_stats

DATA_DIR = os.path.join("data", "soil", "european-soil-database-derived-data")

DATA_FILE = os.path.join(DATA_DIR, "IE_TAWC.tif")

data = rxr.open_rasterio(DATA_FILE, chunks="auto", masked=True)

data

<xarray.DataArray (band: 1, y: 409, x: 418)>
dask.array<open_rasterio-f60863f961dba41a45774b674cfc20e4<this-array>, shape=(1, 409, 418), dtype=float32, chunksize=(1, 409, 418), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 2.918e+06 2.92e+06 2.92e+06 ... 3.334e+06 3.336e+06
  * y            (y) float64 3.728e+06 3.728e+06 ... 3.322e+06 3.32e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    units:          milimeter
    scale_factor:   1.0
    add_offset:     0.0

xarray.DataArray

• band: 1
• y: 409
• x: 418

• dask.array<chunksize=(1, 409, 418), meta=np.ndarray>

  Array Chunk Bytes 667.82 kiB 667.82 kiB Shape (1, 409, 418) (1, 409, 418) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • x
    (x)
    float64
    2.918e+06 2.92e+06 ... 3.336e+06

    array([2918500., 2919500., 2920500., ..., 3333500., 3334500., 3335500.])

  • y
    (y)
    float64
    3.728e+06 3.728e+06 ... 3.32e+06

    array([3728500., 3727500., 3726500., ..., 3322500., 3321500., 3320500.])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    LOCAL_CS["Unknown",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    spatial_ref :

    LOCAL_CS["Unknown",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    GeoTransform :

    2918000.0 1000.0 0.0 3729000.0 0.0 -1000.0

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Int64Index([1], dtype='int64', name='band'))

  • x
    PandasIndex

    PandasIndex(Float64Index([2918500.0, 2919500.0, 2920500.0, 2921500.0, 2922500.0, 2923500.0,
                  2924500.0, 2925500.0, 2926500.0, 2927500.0,
                  ...
                  3326500.0, 3327500.0, 3328500.0, 3329500.0, 3330500.0, 3331500.0,
                  3332500.0, 3333500.0, 3334500.0, 3335500.0],
                 dtype='float64', name='x', length=418))

  • y
    PandasIndex

    PandasIndex(Float64Index([3728500.0, 3727500.0, 3726500.0, 3725500.0, 3724500.0, 3723500.0,
                  3722500.0, 3721500.0, 3720500.0, 3719500.0,
                  ...
                  3329500.0, 3328500.0, 3327500.0, 3326500.0, 3325500.0, 3324500.0,
                  3323500.0, 3322500.0, 3321500.0, 3320500.0],
                 dtype='float64', name='y', length=409))

• Attributes: (4)

  AREA_OR_POINT :

  Area

  units :

  milimeter

  scale_factor :

  1.0

  add_offset :

  0.0

data.rio.crs

CRS.from_wkt('LOCAL_CS["Unknown",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

# ETRS 89 LAEA
data_crs = 3035

data.rio.write_crs(data_crs, inplace=True)

<xarray.DataArray (band: 1, y: 409, x: 418)>
dask.array<open_rasterio-f60863f961dba41a45774b674cfc20e4<this-array>, shape=(1, 409, 418), dtype=float32, chunksize=(1, 409, 418), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 2.918e+06 2.92e+06 2.92e+06 ... 3.334e+06 3.336e+06
  * y            (y) float64 3.728e+06 3.728e+06 ... 3.322e+06 3.32e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    units:          milimeter
    scale_factor:   1.0
    add_offset:     0.0

xarray.DataArray

• band: 1
• y: 409
• x: 418

• dask.array<chunksize=(1, 409, 418), meta=np.ndarray>

  Array Chunk Bytes 667.82 kiB 667.82 kiB Shape (1, 409, 418) (1, 409, 418) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • x
    (x)
    float64
    2.918e+06 2.92e+06 ... 3.336e+06

    array([2918500., 2919500., 2920500., ..., 3333500., 3334500., 3335500.])

  • y
    (y)
    float64
    3.728e+06 3.728e+06 ... 3.32e+06

    array([3728500., 3727500., 3726500., ..., 3322500., 3321500., 3320500.])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314140356

    inverse_flattening :

    298.257222101

    reference_ellipsoid_name :

    GRS 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    ETRS89

    horizontal_datum_name :

    European Terrestrial Reference System 1989

    projected_crs_name :

    ETRS89-extended / LAEA Europe

    grid_mapping_name :

    lambert_azimuthal_equal_area

    latitude_of_projection_origin :

    52.0

    longitude_of_projection_origin :

    10.0

    false_easting :

    4321000.0

    false_northing :

    3210000.0

    spatial_ref :

    PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]

    GeoTransform :

    2918000.0 1000.0 0.0 3729000.0 0.0 -1000.0

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Int64Index([1], dtype='int64', name='band'))

  • x
    PandasIndex

    PandasIndex(Float64Index([2918500.0, 2919500.0, 2920500.0, 2921500.0, 2922500.0, 2923500.0,
                  2924500.0, 2925500.0, 2926500.0, 2927500.0,
                  ...
                  3326500.0, 3327500.0, 3328500.0, 3329500.0, 3330500.0, 3331500.0,
                  3332500.0, 3333500.0, 3334500.0, 3335500.0],
                 dtype='float64', name='x', length=418))

  • y
    PandasIndex

    PandasIndex(Float64Index([3728500.0, 3727500.0, 3726500.0, 3725500.0, 3724500.0, 3723500.0,
                  3722500.0, 3721500.0, 3720500.0, 3719500.0,
                  ...
                  3329500.0, 3328500.0, 3327500.0, 3326500.0, 3325500.0, 3324500.0,
                  3323500.0, 3322500.0, 3321500.0, 3320500.0],
                 dtype='float64', name='y', length=409))

• Attributes: (4)

  AREA_OR_POINT :

  Area

  units :

  milimeter

  scale_factor :

  1.0

  add_offset :

  0.0

data.rio.crs

CRS.from_epsg(3035)

data.rio.resolution()

(1000.0, -1000.0)

# Ireland boundary
GPKG_BOUNDARY = os.path.join("data", "boundaries", "boundaries_all.gpkg")
ie = gpd.read_file(GPKG_BOUNDARY, layer="NUTS_RG_01M_2021_2157_IE")

ie.crs

<Derived Projected CRS: EPSG:2157>
Name: IRENET95 / Irish Transverse Mercator
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Ireland - onshore. United Kingdom (UK) - Northern Ireland (Ulster) - onshore.
- bounds: (-10.56, 51.39, -5.34, 55.43)
Coordinate Operation:
- name: Irish Transverse Mercator
- method: Transverse Mercator
Datum: IRENET95
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

data.max().values

array(117.38613, dtype=float32)

data.min().values

array(0., dtype=float32)

fig = data.plot(
    robust=True,
    cmap="viridis_r",
    figsize=(7, 7),
    levels=10,
    cbar_kwargs={"label": "Total available water content [mm]"},
)
ie.to_crs(data_crs).boundary.plot(
    ax=fig.axes, color="darkslategrey", linewidth=1
)
plt.title(None)
fig.axes.tick_params(labelbottom=False, labelleft=False)
plt.xlabel(None)
plt.ylabel(None)
plt.tight_layout()
plt.axis("equal")
plt.show()

Grid cells

grid_cells = gpd.read_file(
    os.path.join("data", "ModVege", "params.gpkg"), layer="mera"
)

grid_cells.head()

	x	y	sr	ni	geometry
0	415000.0	497500.0	0.917024	0.35	POLYGON ((415227.037 594343.036, 414771.907 59...
1	417500.0	460000.0	0.763591	0.35	POLYGON ((424513.435 557931.477, 424058.108 56...
2	417500.0	462500.0	0.763591	0.35	POLYGON ((424058.108 560389.011, 423602.795 56...
3	417500.0	492500.0	0.917024	0.35	POLYGON ((418595.321 589882.162, 418140.190 59...
4	417500.0	495000.0	0.917024	0.35	POLYGON ((418140.190 592340.144, 417685.076 59...

grid_cells.crs

<Derived Projected CRS: EPSG:2157>
Name: IRENET95 / Irish Transverse Mercator
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Ireland - onshore. United Kingdom (UK) - Northern Ireland (Ulster) - onshore.
- bounds: (-10.56, 51.39, -5.34, 55.43)
Coordinate Operation:
- name: Irish Transverse Mercator
- method: Transverse Mercator
Datum: IRENET95
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

grid_cells.shape

(14490, 5)

fig = data.plot(
    robust=True,
    cmap="viridis_r",
    figsize=(7, 7),
    levels=10,
    cbar_kwargs={"label": "Total available water content [mm]"},
)
grid_cells.to_crs(data_crs).boundary.plot(
    ax=fig.axes, color="darkslategrey", linewidth=1
)
plt.title(None)
fig.axes.tick_params(labelbottom=False, labelleft=False)
plt.xlabel(None)
plt.ylabel(None)
plt.tight_layout()
plt.axis("equal")
plt.show()

Zonal stats

grid_cells = gpd.GeoDataFrame.from_features(
    zonal_stats(
        vectors=grid_cells.to_crs(data_crs),
        raster=os.path.join(DATA_DIR, "IE_TAWC.tif"),
        stats=["count", "mean"],
        geojson_out=True,
        nodata=-999999,
    ),
    crs=data_crs,
).to_crs(grid_cells.crs)

grid_cells.head()

	geometry	x	y	sr	ni	mean	count
0	POLYGON ((415227.037 594343.036, 414771.907 59...	415000.0	497500.0	0.917024	0.35	0.000000	8
1	POLYGON ((424513.435 557931.477, 424058.108 56...	417500.0	460000.0	0.763591	0.35	0.000000	6
2	POLYGON ((424058.108 560389.011, 423602.795 56...	417500.0	462500.0	0.763591	0.35	0.000000	8
3	POLYGON ((418595.321 589882.162, 418140.190 59...	417500.0	492500.0	0.917024	0.35	8.697089	9
4	POLYGON ((418140.190 592340.144, 417685.076 59...	417500.0	495000.0	0.917024	0.35	0.000000	6

grid_cells.shape

(14490, 7)

grid_cells["mean"].min()

0.0

grid_cells["mean"].max()

117.38613552517361

grid_cells["count"].min()

4

grid_cells["count"].max()

9

grid_cells[grid_cells["count"] == 0]

	geometry	x	y	sr	ni	mean	count

axs = grid_cells.plot(
    column="mean",
    cmap="Spectral",
    scheme="equal_interval",
    edgecolor="darkslategrey",
    linewidth=0.2,
    figsize=(6, 7),
    legend=True,
    legend_kwds={"loc": "upper left", "fmt": "{:.2f}", "title": "WHC [mm]"},
    missing_kwds={
        "color": "darkslategrey",
        "edgecolor": "darkslategrey",
        "label": "No data",
    },
)
for legend_handle in axs.get_legend().legend_handles:
    legend_handle.set_markeredgewidth(0.2)
    legend_handle.set_markeredgecolor("darkslategrey")
axs.tick_params(labelbottom=False, labelleft=False)
plt.axis("equal")
plt.tight_layout()
plt.show()

grid_cells["whc"] = grid_cells["mean"]

grid_cells.drop(columns=["mean", "count"], inplace=True)

grid_cells.head()

	geometry	x	y	sr	ni	whc
0	POLYGON ((415227.037 594343.036, 414771.907 59...	415000.0	497500.0	0.917024	0.35	0.000000
1	POLYGON ((424513.435 557931.477, 424058.108 56...	417500.0	460000.0	0.763591	0.35	0.000000
2	POLYGON ((424058.108 560389.011, 423602.795 56...	417500.0	462500.0	0.763591	0.35	0.000000
3	POLYGON ((418595.321 589882.162, 418140.190 59...	417500.0	492500.0	0.917024	0.35	8.697089
4	POLYGON ((418140.190 592340.144, 417685.076 59...	417500.0	495000.0	0.917024	0.35	0.000000

grid_cells.to_file(
    os.path.join("data", "ModVege", "params.gpkg"), layer="mera"
)