Nitrogen nutritional index - gridded for MÉRA

Soil chemical properties based on LUCAS topsoil data (Ballabio et al., 2019; European Commission, n.d.; Panagos et al., 2022; Panagos et al., 2012): https://esdac.jrc.ec.europa.eu/content/chemical-properties-european-scale-based-lucas-topsoil-data

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

DATA_DIR = os.path.join(
    "data",
    "soil",
    "chemical-properties-european-scale-based-lucas-topsoil-data",
)

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

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

data

<xarray.DataArray (band: 1, y: 818, x: 835)>
dask.array<open_rasterio-9a5b296a8d3bf80d52b3a819e76bdfdf<this-array>, shape=(1, 818, 835), dtype=float32, chunksize=(1, 818, 835), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 2.918e+06 2.919e+06 ... 3.335e+06 3.335e+06
  * y            (y) float64 3.728e+06 3.728e+06 3.727e+06 ... 3.32e+06 3.32e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:           Area
    DataType:                Generic
    RepresentationType:      ATHEMATIC
    STATISTICS_COVARIANCES:  0.9927304438511346
    STATISTICS_MAXIMUM:      15.40726184845
    STATISTICS_MEAN:         2.1118209123088
    STATISTICS_MINIMUM:      0
    STATISTICS_SKIPFACTORX:  1
    STATISTICS_SKIPFACTORY:  1
    STATISTICS_STDDEV:       0.99635859199945
    scale_factor:            1.0
    add_offset:              0.0

xarray.DataArray

• band: 1
• y: 818
• x: 835

• dask.array<chunksize=(1, 818, 835), meta=np.ndarray>

  Array Chunk Bytes 2.61 MiB 2.61 MiB Shape (1, 818, 835) (1, 818, 835) 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.919e+06 ... 3.335e+06

    array([2918450., 2918950., 2919450., ..., 3334450., 3334950., 3335450.])

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

    array([3728450., 3727950., 3727450., ..., 3320950., 3320450., 3319950.])

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

    2918199.9999999995 500.0 0.0 3728700.0 0.0 -500.0

    array(0)

• Indexes: (3)
  • band
    PandasIndex

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

  • x
    PandasIndex

    PandasIndex(Float64Index([2918449.9999999995, 2918949.9999999995, 2919449.9999999995,
                  2919949.9999999995, 2920449.9999999995, 2920949.9999999995,
                  2921449.9999999995, 2921949.9999999995, 2922449.9999999995,
                  2922949.9999999995,
                  ...
                  3330949.9999999995, 3331449.9999999995, 3331949.9999999995,
                  3332449.9999999995, 3332949.9999999995, 3333449.9999999995,
                  3333949.9999999995, 3334449.9999999995, 3334949.9999999995,
                  3335449.9999999995],
                 dtype='float64', name='x', length=835))

  • y
    PandasIndex

    PandasIndex(Float64Index([3728450.0, 3727950.0, 3727450.0, 3726950.0, 3726450.0, 3725950.0,
                  3725450.0, 3724950.0, 3724450.0, 3723950.0,
                  ...
                  3324450.0, 3323950.0, 3323450.0, 3322950.0, 3322450.0, 3321950.0,
                  3321450.0, 3320950.0, 3320450.0, 3319950.0],
                 dtype='float64', name='y', length=818))

• Attributes: (12)

  AREA_OR_POINT :

  Area

  DataType :

  Generic

  RepresentationType :

  ATHEMATIC

  STATISTICS_COVARIANCES :

  0.9927304438511346

  STATISTICS_MAXIMUM :

  15.40726184845

  STATISTICS_MEAN :

  2.1118209123088

  STATISTICS_MINIMUM :

  0

  STATISTICS_SKIPFACTORX :

  1

  STATISTICS_SKIPFACTORY :

  1

  STATISTICS_STDDEV :

  0.99635859199945

  scale_factor :

  1.0

  add_offset :

  0.0

data.rio.crs

CRS.from_epsg(3035)

data.rio.resolution()

(500.0, -500.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(9.951199, dtype=float32)

data.min().values

array(1.2719605, dtype=float32)

fig = data.plot(
    robust=True,
    cmap="viridis_r",
    figsize=(7, 7),
    levels=6,
    cbar_kwargs={"label": "Topsoil nitrogen content [g kg⁻¹]"},
)
ie.to_crs(data.rio.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	geometry
0	415000.0	497500.0	0.917024	POLYGON ((415227.037 594343.036, 414771.907 59...
1	417500.0	460000.0	0.763591	POLYGON ((424513.435 557931.478, 424058.108 56...
2	417500.0	462500.0	0.763591	POLYGON ((424058.108 560389.011, 423602.795 56...
3	417500.0	492500.0	0.917024	POLYGON ((418595.321 589882.162, 418140.190 59...
4	417500.0	495000.0	0.917024	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, 4)

fig = data.plot(
    robust=True,
    cmap="viridis_r",
    figsize=(7, 7),
    levels=6,
    cbar_kwargs={"label": "Topsoil nitrogen content [g kg⁻¹]"},
)
grid_cells.to_crs(data.rio.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.rio.crs),
        raster=os.path.join(DATA_DIR, "IE_N.tif"),
        stats=["count", "mean"],
        geojson_out=True,
    ),
    crs=data.rio.crs,
).to_crs(grid_cells.crs)

grid_cells.head()

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

grid_cells.shape

(14490, 6)

grid_cells["mean"].min()

1.7392230224609375

grid_cells["mean"].max()

8.040995597839355

grid_cells["count"].min()

0

grid_cells["count"].max()

26

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

	geometry	x	y	sr	count	mean
0	POLYGON ((415227.037 594343.036, 414771.907 59...	415000.0	497500.0	0.917024	0	NaN
1	POLYGON ((424513.435 557931.477, 424058.108 56...	417500.0	460000.0	0.763591	0	NaN
2	POLYGON ((424058.108 560389.011, 423602.795 56...	417500.0	462500.0	0.763591	0	NaN
3	POLYGON ((418595.321 589882.162, 418140.190 59...	417500.0	492500.0	0.917024	0	NaN
4	POLYGON ((418140.190 592340.144, 417685.076 59...	417500.0	495000.0	0.917024	0	NaN
...	...	...	...	...	...	...
14483	POLYGON ((759655.600 876544.928, 759200.942 87...	805000.0	712500.0	1.581323	0	NaN
14484	POLYGON ((759200.942 879003.081, 758746.267 88...	805000.0	715000.0	0.900000	0	NaN
14485	POLYGON ((758746.267 881461.207, 758291.577 88...	805000.0	717500.0	2.021384	0	NaN
14486	POLYGON ((758291.577 883919.305, 757836.870 88...	805000.0	720000.0	2.021384	0	NaN
14487	POLYGON ((765750.631 857333.286, 765296.086 85...	807500.0	692500.0	2.497536	0	NaN

729 rows × 6 columns

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": "Topsoil N [g kg⁻¹]",
    },
    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()

Normalise

# normalise between 1.0 and 0.35
grid_cells["ni"] = 0.35 + (
    (grid_cells["mean"] - float(grid_cells["mean"].min())) * (1.0 - 0.35)
) / (float(grid_cells["mean"].max()) - float(grid_cells["mean"].min()))

grid_cells.head()

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

grid_cells["ni"].max()

1.0

grid_cells["ni"].min()

0.35

axs = grid_cells.plot(
    column="ni",
    cmap="Spectral_r",
    scheme="equal_interval",
    edgecolor="darkslategrey",
    linewidth=0.2,
    figsize=(6, 7),
    legend=True,
    legend_kwds={"loc": "upper left", "fmt": "{:.2f}", "title": "NNI"},
    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["ni"] = grid_cells["ni"].fillna(grid_cells["ni"].min())

grid_cells.head()

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

axs = grid_cells.plot(
    column="ni",
    cmap="Spectral_r",
    scheme="equal_interval",
    edgecolor="darkslategrey",
    linewidth=0.2,
    figsize=(6, 7),
    legend=True,
    legend_kwds={"loc": "upper left", "fmt": "{:.2f}", "title": "NNI"},
    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.to_file(
    os.path.join("data", "ModVege", "params.gpkg"), layer="mera"
)