Merge branch 'aggregation_plotter' into 'master'

added mnp_hokaggregation.py

See merge request !6

benb.py

@@ -27,6 +27,66 @@ def rand_web_color_hex():
     return "#{}".format(rgb)
 
 
+def fix_bt(code_in, as_mnp=False):
+    """
+    Parser for Beheertype codes to repair abbreviated codes and/or MNP extension
+    eg. 02.01     --> N02.01 if as_mnp==False
+        02.01     --> N01.01.00 if as_mnp==True
+        N05.01    --> N05.01  if as_mnp==False
+        N05.01    --> N05.01.00  if as_mnp==True
+        N04.02.01 --> N04.02.01 if as_mnp==False
+        N04.02.01 --> N04.02.01 if as_mnp==True
+    BT Code is as follows <N><AA>.<BB>.<CC>
+    N = letter, may be missing
+    A = two digits indicating top level beheertype "TOP"
+    B = two digits indicating sub level beheertype "SUB"
+    C = two digits indicating neergeschaald beheertype if > 00. Else "00" as used in MNP param notation  "NEER"
+    N04.02.00 is equivalent to N04.02
+    :param code_in: beheertypecode in any format
+    :param as_mnp: return MNP style beheertypecode, default=False
+    :return: beheertypecode
+    """
+    if isinstance(code_in, numbers.Number):
+        code_in = str(code_in)
+
+    # Identify parts of code
+    parts = code_in.split(".")
+    if len(parts) == 1:
+        raise Exception("Unexpected BT code: {}".format(code_in))
+    elif len(parts) == 2:
+        top, sub = parts
+        neer = None
+    elif len(parts) == 3:
+        top, sub, neer = parts
+
+    # Verify TOP and add "N" is required
+    if top.isdigit():
+        top = "N{0}".format(top.zfill(2))
+    else:
+        assert top[0].isupper(), "Unexpected BT Code format: {}".format(code_in)
+
+    # Verify SUB
+    assert sub.isdigit and len(sub) == 2, "Unexpected BT Code format: {}".format(
+        code_in
+    )
+
+    # Verify NEER and check if it is neergeschaald or not
+    keep_neer = False
+    if neer:
+        assert neer.isdigit and len(neer) == 2, "Unexpected BT Code format: {}".format(
+            code_in
+        )
+        if int(neer) > 0:
+            keep_neer = True
+
+    # Construct output
+    head = "{}.{}".format(top, sub)
+    if keep_neer:
+        return "{}.{}".format(head, neer)
+    elif as_mnp:
+        return "{}.00".format(head)
+    else:
+        return head
 
 
 def pnt2square_id(east_north, size, of="full"):
@@ -69,6 +129,21 @@ def pnt2square_id(east_north, size, of="full"):
         raise
 
 
+def draagkracht_ww(frequentie, trouw, th):
+    """
+    Draagkracht formule volgens Wieger Wamelink. Zie MS Teams discussie 16-07-2020 MNP2020/Fase 7 Objectivering Draagk..
+    :param frequentie: gemiddelde frequentie van soort X binnen beheertype Y
+    :param trouw: gemiddelde trouw van soort X binnen beheertype Y
+    :param th: minumum draagkracht value
+    :return: draagkracht soort X - beheertype Y combinatie, gemaximaliseerd tot 1.
+    """
+    dk = np.min([np.multiply(np.divide(np.max([frequentie, trouw]), 100), 5), 1])
+    if dk >= th:
+        return dk
+    else:
+        return np.nan
+
+
 def gen_squares(x_ul, y_ul, nrow, ncol, size):
     """
     Generate geodataframe with square polygons

mnp_hokaggregation.py

+import rasterio as rio
+from rasterio.plot import show
+import pandas as pd
+import geopandas as gp
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import os
+
+
+def mnp_hokaggregation(hok_csv, hokken_polygon_path, background_raster_path, outputdir):
+    """
+    make aggregationplots for use in manual validation of mnp runs. can both aggregate validation data and species
+    model gridcells per hok.
+    :param hok_csv: csv file with value per hok
+    :param hokken_polygon_path: path to hokken polygon file
+    :param background_raster_path: path to raster file to be used as background. raster should only contain on value
+     lower than or equal to 0, and 1's.
+    :param outputdir: directory to save plot
+    :return:
+    """
+
+    # load hokken
+    hokken = gp.read_file(hokken_polygon_path)
+    if "ID" not in hokken.columns:
+        raise KeyError("No hok identifiercolumn called ID column in shapefile")
+
+    # load and prepare background raster
+    rast = rio.open(background_raster_path)
+    background = rast.read(1)
+    background[background < 0] = 0
+    cmp = colors.ListedColormap([[0.97, 0.97, 0.97, 1], [0.8, 0.8, 0.8, 1]])
+
+    # Load counts data and shapes
+    counts = pd.read_csv(hok_csv)
+    colnames = list(counts.columns)
+    colnames[colnames.index("UH_ID")] = "ID"
+    counts.columns = colnames
+    counthokken = pd.merge(hokken, counts, on="ID")
+
+    if "COUNT_KEY_VALUE" in colnames:
+        bounds = np.array([1, 10, 100, 1000, 10000, 40000])
+        explabel = "_"
+        valuecolumn = "COUNT_KEY_VALUE"
+        extend = "neither"
+        ylabel_text = "Number of gridcells in Hok"
+
+    elif "OPS_COUNT" in colnames:
+        bounds = np.array([1, 10, 50, 100, 500, 1000])
+        explabel = "_validation"
+        valuecolumn = "OPS_COUNT"
+        extend = "max"
+        ylabel_text = "number of observations in NDFF"
+
+    else:
+        raise Exception("no value column in columns")
+
+    # Make the figure and plot, you know the drill
+    plt.ioff()
+    fig, ax = plt.subplots(figsize=(15, 15))
+    plt.axis("off")
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes("right", size=0.25, pad=0.1)
+    show(background, transform=rast.transform, ax=ax, cmap=cmp)
+    counthokken.plot(
+        ax=ax,
+        alpha=0.8,
+        column=valuecolumn,
+        cmap="YlOrRd",
+        edgecolor="darkgray",
+        norm=colors.BoundaryNorm(boundaries=bounds, ncolors=256, extend=extend),
+        linewidth=0.8,
+        legend=True,
+        cax=cax,
+    )
+    # cax is a seperate axes dedicated tot the colorbar legend on the right
+    cax.set_ylabel(ylabel_text, size=20, fontweight="semibold")
+    cax.tick_params(labelsize=15)
+
+    # make outputdir if not exists and save figure
+    outputdir = outputdir
+    if not os.path.exists(outputdir):
+        os.makedirs(outputdir)
+    plt.savefig(
+        outputdir + "/" + os.path.basename(hok_csv)[:-4] + ".png", bbox_inches="tight", dpi=400
+    )
+    fig.clear()
+    plt.close(fig)