Soil water-holding capacity - gridded for EURO-CORDEX

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
from zipfile import BadZipFile, ZipFile
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")

ZIP_FILE = os.path.join(DATA_DIR, "STU_EU_Layers.zip")

# list of files/folders in the ZIP archive
ZipFile(ZIP_FILE).namelist()

['SMU_EU_S_TAWC.RDC',
 'SMU_EU_S_TAWC.rst',
 'SMU_EU_T_TAWC.RDC',
 'SMU_EU_T_TAWC.rst',
 'STU_EU_ALLOCATION.RDC',
 'STU_EU_ALLOCATION.rst',
 'STU_EU_DEPTH_ROOTS.RDC',
 'STU_EU_DEPTH_ROOTS.rst',
 'STU_EU_S_BD.RDC',
 'STU_EU_S_BD.rst',
 'STU_EU_S_CLAY.RDC',
 'STU_EU_S_CLAY.rst',
 'STU_EU_S_GRAVEL.RDC',
 'STU_EU_S_GRAVEL.rst',
 'STU_EU_S_OC.RDC',
 'STU_EU_S_OC.rst',
 'STU_EU_S_SAND.RDC',
 'STU_EU_S_SAND.rst',
 'STU_EU_S_SILT.RDC',
 'STU_EU_S_SILT.rst',
 'STU_EU_S_TAWC.RDC',
 'STU_EU_S_TAWC.rst',
 'STU_EU_S_TEXT_CLS.RDC',
 'STU_EU_S_TEXT_CLS.rst',
 'STU_EU_T_BD.RDC',
 'STU_EU_T_BD.rst',
 'STU_EU_T_CLAY.RDC',
 'STU_EU_T_CLAY.rst',
 'STU_EU_T_GRAVEL.RDC',
 'STU_EU_T_GRAVEL.rst',
 'STU_EU_T_OC.RDC',
 'STU_EU_T_OC.rst',
 'STU_EU_T_SAND.RDC',
 'STU_EU_T_SAND.rst',
 'STU_EU_T_SILT.RDC',
 'STU_EU_T_SILT.rst',
 'STU_EU_T_TAWC.RDC',
 'STU_EU_T_TAWC.rst',
 'STU_EU_T_TEXT_CLS.RDC',
 'STU_EU_T_TEXT_CLS.rst']

# extract the archive
try:
    z = ZipFile(ZIP_FILE)
    z.extractall(DATA_DIR)
except BadZipFile:
    print("There were issues with the file", ZIP_FILE)

DATA_FILE = os.path.join(DATA_DIR, "STU_EU_T_TAWC.rst")

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

data

<xarray.DataArray (band: 1, y: 4600, x: 5900)>
dask.array<open_rasterio-fe4b1480115aa1d61925a510437a7116<this-array>, shape=(1, 4600, 5900), dtype=float32, chunksize=(1, 4600, 5900), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 1.5e+06 1.502e+06 1.502e+06 ... 7.398e+06 7.4e+06
  * y            (y) float64 5.5e+06 5.498e+06 5.498e+06 ... 9.015e+05 9.005e+05
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    units:         milimeter

xarray.DataArray

• band: 1
• y: 4600
• x: 5900

• dask.array<chunksize=(1, 4600, 5900), meta=np.ndarray>

  Array Chunk Bytes 103.53 MiB 103.53 MiB Shape (1, 4600, 5900) (1, 4600, 5900) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

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

    array([1])

  • x
    (x)
    float64
    1.5e+06 1.502e+06 ... 7.4e+06

    array([1500500., 1501500., 1502500., ..., 7397500., 7398500., 7399500.])

  • y
    (y)
    float64
    5.5e+06 5.498e+06 ... 9.005e+05

    array([5499500., 5498500., 5497500., ...,  902500.,  901500.,  900500.])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    LOCAL_CS["Unknown",UNIT["Meter",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    spatial_ref :

    LOCAL_CS["Unknown",UNIT["Meter",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    GeoTransform :

    1500000.0 1000.0 0.0 5500000.0 0.0 -1000.0

    array(0)

• Indexes: (3)
  • band
    PandasIndex

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

  • x
    PandasIndex

    PandasIndex(Float64Index([1500500.0, 1501500.0, 1502500.0, 1503500.0, 1504500.0, 1505500.0,
                  1506500.0, 1507500.0, 1508500.0, 1509500.0,
                  ...
                  7390500.0, 7391500.0, 7392500.0, 7393500.0, 7394500.0, 7395500.0,
                  7396500.0, 7397500.0, 7398500.0, 7399500.0],
                 dtype='float64', name='x', length=5900))

  • y
    PandasIndex

    PandasIndex(Float64Index([5499500.0, 5498500.0, 5497500.0, 5496500.0, 5495500.0, 5494500.0,
                  5493500.0, 5492500.0, 5491500.0, 5490500.0,
                  ...
                   909500.0,  908500.0,  907500.0,  906500.0,  905500.0,  904500.0,
                   903500.0,  902500.0,  901500.0,  900500.0],
                 dtype='float64', name='y', length=4600))

• Attributes: (3)

  scale_factor :

  1.0

  add_offset :

  0.0

  units :

  milimeter

data.rio.crs

CRS.from_wkt('LOCAL_CS["Unknown",UNIT["Meter",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

# manually assign CRS: ETRS 89 LAEA
data_crs = 3035

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

<xarray.DataArray (band: 1, y: 4600, x: 5900)>
dask.array<open_rasterio-fe4b1480115aa1d61925a510437a7116<this-array>, shape=(1, 4600, 5900), dtype=float32, chunksize=(1, 4600, 5900), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 1.5e+06 1.502e+06 1.502e+06 ... 7.398e+06 7.4e+06
  * y            (y) float64 5.5e+06 5.498e+06 5.498e+06 ... 9.015e+05 9.005e+05
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
    units:         milimeter

xarray.DataArray

• band: 1
• y: 4600
• x: 5900

• dask.array<chunksize=(1, 4600, 5900), meta=np.ndarray>

  Array Chunk Bytes 103.53 MiB 103.53 MiB Shape (1, 4600, 5900) (1, 4600, 5900) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

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

    array([1])

  • x
    (x)
    float64
    1.5e+06 1.502e+06 ... 7.4e+06

    array([1500500., 1501500., 1502500., ..., 7397500., 7398500., 7399500.])

  • y
    (y)
    float64
    5.5e+06 5.498e+06 ... 9.005e+05

    array([5499500., 5498500., 5497500., ...,  902500.,  901500.,  900500.])

  • 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 :

    1500000.0 1000.0 0.0 5500000.0 0.0 -1000.0

    array(0)

• Indexes: (3)
  • band
    PandasIndex

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

  • x
    PandasIndex

    PandasIndex(Float64Index([1500500.0, 1501500.0, 1502500.0, 1503500.0, 1504500.0, 1505500.0,
                  1506500.0, 1507500.0, 1508500.0, 1509500.0,
                  ...
                  7390500.0, 7391500.0, 7392500.0, 7393500.0, 7394500.0, 7395500.0,
                  7396500.0, 7397500.0, 7398500.0, 7399500.0],
                 dtype='float64', name='x', length=5900))

  • y
    PandasIndex

    PandasIndex(Float64Index([5499500.0, 5498500.0, 5497500.0, 5496500.0, 5495500.0, 5494500.0,
                  5493500.0, 5492500.0, 5491500.0, 5490500.0,
                  ...
                   909500.0,  908500.0,  907500.0,  906500.0,  905500.0,  904500.0,
                   903500.0,  902500.0,  901500.0,  900500.0],
                 dtype='float64', name='y', length=4600))

• Attributes: (3)

  scale_factor :

  1.0

  add_offset :

  0.0

  units :

  milimeter

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

# clip raster to Ireland's boundary
data = rxr.open_rasterio(DATA_FILE, chunks="auto", masked=True).rio.clip(
    ie.buffer(6500).to_crs(data_crs)
)

data

<xarray.DataArray (band: 1, y: 409, x: 418)>
dask.array<getitem, 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:
    scale_factor:  1.0
    add_offset:    0.0
    units:         milimeter

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 5 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

    axis :

    X

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

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

    axis :

    Y

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

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    LOCAL_CS["Unknown",UNIT["Meter",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    spatial_ref :

    LOCAL_CS["Unknown",UNIT["Meter",1],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: (3)

  scale_factor :

  1.0

  add_offset :

  0.0

  units :

  milimeter

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("")
plt.ylabel("")
plt.tight_layout()
plt.axis("equal")
plt.show()

# export to GeoTIFF
data.rio.to_raster(os.path.join(DATA_DIR, "IE_TAWC.tif"))

Grid cells

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

grid_cells.head()

	rlon	rlat	sr	ni	geometry
0	-17.265	4.345	0.763591	0.35	POLYGON ((419870.326 550490.053, 415715.725 56...
1	-17.265	4.455	0.000000	0.35	POLYGON ((415715.725 562009.784, 411561.600 57...
2	-17.265	4.565	0.000000	0.35	POLYGON ((411561.600 573530.159, 407407.964 58...
3	-17.265	4.675	0.000000	0.35	POLYGON ((407407.964 585051.187, 403254.829 59...
4	-17.265	4.785	0.000000	0.35	POLYGON ((403254.829 596572.878, 399102.207 60...

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

(765, 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("")
plt.ylabel("")
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	rlon	rlat	sr	ni	mean	count
0	POLYGON ((419870.326 550490.053, 415715.725 56...	-17.265	4.345	0.763591	0.35	0.0	119
1	POLYGON ((415715.725 562009.784, 411561.600 57...	-17.265	4.455	0.000000	0.35	0.0	45
2	POLYGON ((411561.600 573530.159, 407407.964 58...	-17.265	4.565	0.000000	0.35	0.0	54
3	POLYGON ((407407.964 585051.187, 403254.829 59...	-17.265	4.675	0.000000	0.35	0.0	116
4	POLYGON ((403254.829 596572.878, 399102.207 60...	-17.265	4.785	0.000000	0.35	0.0	67

grid_cells.shape

(765, 7)

grid_cells["mean"].min()

0.0

grid_cells["mean"].max()

114.76114198825503

grid_cells["count"].min()

0

grid_cells["count"].max()

154

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

	geometry	rlon	rlat	sr	ni	mean	count
99	POLYGON ((474491.495 687204.081, 470328.825 69...	-16.385	5.335	0.0	0.35	NaN	0
458	POLYGON ((556849.412 925216.933, 552662.092 93...	-14.955	7.095	0.0	0.35	NaN	0
611	POLYGON ((734722.641 651078.195, 730526.414 66...	-14.295	4.235	0.0	0.35	NaN	0
731	POLYGON ((753193.890 814183.944, 748979.265 82...	-13.635	5.555	0.0	0.35	NaN	0

axs = grid_cells.plot(
    column="mean",
    cmap="Spectral_r",
    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()

# fill no data with min value
grid_cells["mean"] = grid_cells["mean"].fillna(grid_cells["mean"].min())

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

axs = grid_cells.plot(
    column="mean",
    cmap="Spectral_r",
    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.drop(columns=["mean", "count"], inplace=True)

grid_cells.head()

	geometry	rlon	rlat	sr	ni	whc
0	POLYGON ((419870.326 550490.053, 415715.725 56...	-17.265	4.345	0.763591	0.35	0.0
1	POLYGON ((415715.725 562009.784, 411561.600 57...	-17.265	4.455	0.000000	0.35	0.0
2	POLYGON ((411561.600 573530.159, 407407.964 58...	-17.265	4.565	0.000000	0.35	0.0
3	POLYGON ((407407.964 585051.187, 403254.829 59...	-17.265	4.675	0.000000	0.35	0.0
4	POLYGON ((403254.829 596572.878, 399102.207 60...	-17.265	4.785	0.000000	0.35	0.0

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