Committing from desktop new DK scripts

sample/dk.py

 """
-Determine draagkracht for species
+Determine draagkracht for vogel and vlinder using NDFF database
 Hans Roelofsen, mei 2020
 """
 
@@ -11,16 +11,27 @@ import rasterstats as rast
 import rasterio as rio
 from mnp import helpers as hlp
 # from sample.mnp import helpers as hlp
+# sp = hlp.Species('Tapuit')
+# methode = 'centroid'
 
-# For which species?
 parser = argparse.ArgumentParser()
 parser.add_argument('species', type=str, help='species name')
+parser.add_argument('method', choices=['orig', 'norm', 'centroid'])
 parser.add_argument('--totals', action='store_true', default=False, help='include total hectares of beheertypen?')
 parser.add_argument('--indx', action='store_true', default=False, help='include index?')
+parser.add_argument('--mnp', action='store_true', default=False, help='include MNP draagkracht referentie?')
+parser.add_argument('--pxls', help='output in pixel count instead of hectares', action='store_true', default=False)
+parser.add_argument('--ref_plys', help='path to referentie polygonen')
+
+# Parse arguments
 args = parser.parse_args()
 sp = hlp.Species(args.species)
 totals = args.totals
+methode = args.method
+output_factor = 1 if args.pxls else 0.0625
+output_unit = 'pixels' if args.pxls else 'hectares'
 indx = args.indx
+mnp = args.mnp
 
 # Path to beheertypenkaart
 bt_tiff = r'c:\apps\temp_geodata\nbt\nbt_20190730_111603.tif'
@@ -28,6 +39,8 @@ bt_kaart = rio.open(bt_tiff, mode='r+')
 
 # Table with pixel values, beheertypecodes etc
 nbts = hlp.nbts()
+pxl_2_btB = dict(zip(nbts.pxl, nbts.bt_B))
+btA_2_btB = dict(zip(nbts.bt_A, nbts.bt_B))
 
 # Read NDFF shapefile
 ndff = hlp.read_ndff_shp(sp.ndff_src)
@@ -40,46 +53,69 @@ sel = set.intersection(*[set(queries[x]) for x in queries.keys()])  # for query
 assert len(sel) > 0, 'no observations remaining'
 print('Found {0} NDFF observations for {1}'.format(len(sel), sp.name))
 
-# normalize NDFF observations to 10m radius circle around centroid
-ndff_norm = gp.GeoDataFrame({'ID': [1]*len(sel), 'soort': [sp.name]*len(sel),
-                            'geometry': ndff.loc[sel].centroid.buffer(12.5)}).dissolve(by='ID')
+# Create GDF from NDFF DF following requested method
+if methode == 'norm':
+    ndff_gdf = gp.GeoDataFrame({'ID': [1] * len(sel), 'soort': [sp.name] * len(sel),
+                                'geometry': ndff.loc[sel].centroid.buffer(12.5)}).dissolve(by='ID')
+elif methode == 'centroid':
+    ndff_gdf = gp.GeoDataFrame({'ID': [1] * len(sel), 'soort': [sp.name] * len(sel),
+                                'geometry': ndff.loc[sel].centroid})
+else:  # methode == orig
+    ndff_gdf = gp.GeoDataFrame({'ID': [1] * len(sel), 'soort': [sp.name] * len(sel),
+                                'geometry': ndff.loc[sel]}).dissolve(by='ID')
+
+if methode in ['norm', 'orig']:
+    # Overlay normalized NDFF observations with nbt kaart
+    sp_stats = rast.zonal_stats(vectors=ndff_gdf, raster=bt_kaart.read(1), categorical=True, all_touched=False,
+                                category_map=pxl_2_btB, affine=bt_kaart.affine)
+else:  # methode == 'centroids':
+    # Couple each NDFF observation to a Col-Row index from the bt_kaart using the Affine transformation.
+    # Then drop duplicates to ensure only 1 NDFF observation is counted per bt_kaart pixel.
+    # See: https://www.perrygeo.com/python-affine-transforms.html
 
-# Overlay normalized NDFF observations with nbt kaart
-cat_map = dict(zip(nbts.pxl, nbts.bt_B))
-sp_stats = rast.zonal_stats(vectors=ndff_norm, raster=bt_kaart.read(1), categorical=True, all_touched=False,
-                            category_map=cat_map, affine=bt_kaart.affine)
+    ndff_gdf.loc[:, 'nbt_row_col'] = ndff_gdf.geometry.apply(hlp.xy_2_colrow, rio_object=bt_kaart)
+    ndff_gdf.drop_duplicates(subset='nbt_row_col', inplace=True)
+    sp_stats = rast.zonal_stats(vectors=ndff_gdf, raster=bt_kaart.read(1), categorical=True, all_touched=False,
+                                category_map=pxl_2_btB, affine=bt_kaart.affine)
 
 # Presence of species as pixel count for all beheertypen in de beheertypekaart, dus niet enkel de gevonden.
 species_stats = dict(zip(set(nbts.bt_B), [0] * len(set(nbts.bt_B))))
-for polygon in sp_stats:
-    bts = [s for s in polygon.keys()]
+for waarneming in sp_stats:
+    bts = [s for s in waarneming.keys()]
     for beheertype in bts:
-        species_stats[beheertype] += np.multiply(polygon[beheertype], 0.0625)
+        species_stats[beheertype] += np.multiply(waarneming[beheertype], output_factor)
 
-# Total acreage of each beheertype, assuming pixel size 25m
-bt_u, bt_c = np.unique(bt_kaart.read(1), return_counts=True)
-bt_stats_prelim = dict(zip(bt_u, [np.multiply(bt_count, 0.0625) for bt_count in bt_c]))
-# Note how >1 pixel values can point to 1 beheertype!
-bt_stats = dict(zip(set(nbts.bt_B), np.zeros(nbts.shape[0])))
-for k, v in bt_stats_prelim.items():
-    bt = cat_map[k]
-    bt_stats[bt] += v
-# Check totals of beheertypen!
-assert np.sum([v for _, v in bt_stats.items()]) == np.multiply(np.product(bt_kaart.shape), 0.0625)
+# Total acreage of each beheertype
+if totals:
+    bt_u, bt_c = np.unique(bt_kaart.read(1), return_counts=True)
+    bt_stats_prelim = dict(zip(bt_u, [np.multiply(bt_count, output_factor) for bt_count in bt_c]))
+    # Note how >1 pixel values can point to 1 beheertype!
+    bt_stats = dict(zip(set(nbts.bt_B), np.zeros(nbts.shape[0])))
+    for k, v in bt_stats_prelim.items():
+        bt = pxl_2_btB[k]
+        bt_stats[bt] += v
+    # Check totals of beheertypen!
+    assert np.sum([v for _, v in bt_stats.items()]) == np.multiply(np.product(bt_kaart.shape), output_factor)
+
+# Read draagkrachten as they currently are from MNP parameterset
+if mnp:
+    mnp_stats_prelim = hlp.get_mnp_dk(sp)  # MNP draagkrachten are provided for beheertype_A, ie codes from 20190730
+    mnp_stats = dict(zip(set(nbts.bt_B), np.zeros(nbts.shape[0])))
+    for k, v in mnp_stats_prelim.items():
+        try:
+            mnp_stats[k] += v
+        except KeyError:
+            pass
 
 # Summary dataframe w. Beheertypen as columns and species acreage as rows
+datas = [species_stats]
+indexes = ['{0}_{1}'.format(sp.name, output_unit)]
+if mnp:
+    datas.append(mnp_stats)
+    indexes.append('03_MNP_{0}'.format(sp.name))
 if totals:
-    dk_df = pd.DataFrame([species_stats, bt_stats], index=['{}_ha'.format(sp.name), 'total_ha'])
-else:
-    dk_df = pd.DataFrame(species_stats, index=['{}_ha'.format(sp.name)])
-dk_df.T.to_clipboard(sep=';', index=True if indx else False)
+    datas.append(bt_stats)
+    indexes.append('total_{0}'.format(output_unit))
+pd.DataFrame(data=datas, index=indexes).T.to_clipboard(sep=';', index=True if indx else False)
 print('\n  results are on ze clipboard')
 
-"""
-# Plot bardiagram with BT codes as x-axis labels and areas as % as barheights
-labs = sorted([k for k,_ in species_stats.items()])
-vals = [np.around(np.multiply(np.divide(species_stats[x], np.sum([v for _, v in species_stats.items()])), 100),
-                  decimals=2) for x in labs]
-hlp.plt_bar(labels=labs, values=vals, species=sp.name,
-            ha=np.round(np.multiply(np.sum([v for _, v in species_stats.items()]), 0.0625), 0))
-"""
\ No newline at end of file

sample/dk_adv.py

+"""
+Program for calculating draagkracht (dk) for MNP parameterisation. Sourced from NDFF species observations dataset,
+combined with nominal raster beheertypen (bt). DK is defined between a bt <-> species as the ratio b'teen overal area
+ of a bt and fraction of bt with 1 or more observations of the species.
+
+Variables:
+species_name: name of the species, either vogel or vlinder
+ply_src: polygon shapefile for spatial differentiation of dk (areas)
+nbt_src: path to neergeschaalde beheertypenkaart (bt) (nominal raster)
+
+Method:
+1. total pixels per bt per area
+2. reduce NDFF observations to 1 observation per nbt pixel using forward affine transform
+3. couple each NDFF observation to 1 bt
+4. couple each NDFF observation to 1 area
+5. rework tables to output
+
+                             <SPECIES NAME>
+AreaName|BT1_tot|BT1_presence|BT2_tot|BT2_presence|...|BTn_tot|BTn_presence
+
+By: Hans Roelofsen, WEnR, 29 June 2020
+"""
+
+import argparse
+import os
+import pandas as pd
+import geopandas as gp
+import rasterio as rio
+import rasterstats as rast
+import numpy as np
+from mnp import helpers as hlp
+# from sample.mnp import helpers as hlp
+
+parser = argparse.ArgumentParser()
+parser.add_argument('species', type=str, help='species name')
+
+# Parse arguments
+args = parser.parse_args()
+sp = hlp.Species(args.species)
+
+# Path to beheertypenkaart
+bt_tiff = r'c:\apps\temp_geodata\nbt\nbt_20190730_111603.tif'
+bt_kaart = rio.open(bt_tiff, mode='r+')
+
+# Table with pixel values, beheertypecodes etc
+nbts = hlp.nbts()
+pxl_2_btB = dict(zip(nbts.pxl, nbts.bt_B))
+btA_2_btB = dict(zip(nbts.bt_A, nbts.bt_B))
+
+# Read NDFF shapefile
+ndff = hlp.read_ndff_shp(sp.ndff_src)
+
+# Read areas shapefile
+areas = gp.read_file()
+
+# define queries that isolate NDFF data of interest
+queries = {'q1': ndff.loc[(ndff.year >= 2000) & (ndff.year.notna())].index,
+           'q2': ndff.loc[ndff.nl_name.isin([sp.name])].index,
+           'q3': ndff.loc[ndff.area < 10000].index}
+sel = set.intersection(*[set(queries[x]) for x in queries.keys()])  # for query in selected_queries
+assert len(sel) > 0, 'no observations remaining'
+print('Found {0} NDFF observations for {1}'.format(len(sel), sp.name))
+
+ndff_gdf = gp.GeoDataFrame({'ID': [1] * len(sel), 'soort': [sp.name] * len(sel), 'geometry': ndff.loc[sel].centroid})
+ndff_gdf.loc[:, 'nbt_row_col'] = ndff_gdf.geometry.apply(hlp.xy_2_colrow, rio_object=bt_kaart)
+ndff_gdf.drop_duplicates(subset='nbt_row_col', inplace=True)
+sp_stats = rast.zonal_stats(vectors=ndff_gdf, raster=bt_kaart.read(1), categorical=True, all_touched=False,
+                            category_map=pxl_2_btB, affine=bt_kaart.affine)
+

sample/dk_plant.py

+"""
+Draagkracht calculations for plant species based on DVN vegetatieopnames.
+
+Target format
+<species_code>,<land_type_code>,<land_type_quality>
+
+Species codes and set are sourced from MNP parameter files.
+Relation between species presence and vegetation types is sourced from DVN data, w thanks to Stephan Hennekes
+
+Hans Roelofsen, 23 June 2020
+"""
+
+import os
+import pandas as pd
+
+"""Source data"""
+src_mnp_sp_tax = r'C:\apps\temp_geodata\MNP_DK\07_MNP_versie4_par_population_factors.csv'
+src_mnp_sp_nms = r'c:\apps\temp_geodata\MNP_DK\09_MNP_versie4_Group_Species_valid model_468.csv'
+src_asso_x_bts = r'C:\apps\temp_geodata\MNP_DK\vertaling_IPO_ASSO_2009.txt'
+src_asso_x_sps = r'C:\apps\temp_geodata\MNP_DK\q_syntaxa_soorten_frequentie_trouwsoorten.txt'
+
+mnp_sp_tax = pd.read_csv(src_mnp_sp_tax, sep=',')
+mnp_sp_nms = pd.read_csv(src_mnp_sp_nms, sep=',')
+asso_x_bts = pd.read_csv(src_asso_x_bts, sep=';', comment='#')
+asso_x_sps = pd.read_csv(src_asso_x_sps, sep=';', comment='#', usecols=["Syntaxon_code", "Syntaxon_naam", "Soortnummer",
+                                                                        "Soortnaam", "Frequentie", "Trouw"])
+
+mnp
\ No newline at end of file

sample/mnp/helpers.py

@@ -68,6 +68,7 @@ class Species:
         self.name = sp_name.lower()
         self.soortgr = None
         self.soortnr = None
+        self.mnp_species_code = None
         self.ndff_src = None
         self.nw = None
         self.lpi = None
@@ -90,7 +91,54 @@ class Species:
             self.soortgr = sp_info.loc[sp_info.sp_nm.str.lower() == self.name, 'tax_groep'].item()
             self.soortnr = sp_info.loc[sp_info.sp_nm.str.lower() == self.name, 'sp_nr'].item()
             self.ndff_src = shp_src[self.soortgr]
+            self.mnp_species_code = '{0}{1}'.format({'broedvogel': 'S02', 'vaatplant': 'S09',
+                                                     'dagvlinder': 'S06'}[self.soortgr], str(self.soortnr)[1:])
         except (NameError, ValueError) as e:
             print('{}: {} is not a recognized name'.format(e, self.name))
             raise
 
+
+def colrow_2_xy(colrow, rio_object):
+    """
+    https://www.perrygeo.com/python-affine-transforms.html
+    :param rio_object:
+    :param col: column index
+    :param row: row index
+    :return: easting norhting in CRS
+    """
+    try:
+        a = rio_object.affine
+        col, row = colrow
+        x, y = a * (col, row)
+        return x, y
+    except AttributeError:
+        print("{0} is not a valid rasterio object".format(rio_object))
+        raise
+
+
+def xy_2_colrow(point, rio_object):
+    """
+    https://www.perrygeo.com/python-affine-transforms.html
+    """
+    try:
+        a = rio_object.affine
+        x, y = point.x, point.y
+        col, row = ~a * (x, y)
+        return '{0}_{1}'.format(int(col//1), int(row//1))
+    except AttributeError as e:
+        print(e)
+        raise
+        # is dit nuttig?
+
+
+def get_mnp_dk(species):
+    try:
+        df = pd.read_csv(r'c:\apps\temp_geodata\MNP_DK\03_MNP_versie6_par_density_factors_BT2019_v2.csv', sep=',',
+                         usecols=['Species_code', 'Land_type_code', 'Land_type_quality'])
+        df['bt_B'] = np.where(df.Land_type_code.str.slice(7) == '00', df.Land_type_code.str.slice(0, 6),
+                              df.Land_type_code)
+        foo = df.loc[df.Species_code == species.mnp_species_code, ['bt_B', 'Land_type_quality']]
+        return dict(zip(foo.bt_B, foo.Land_type_quality))
+    except AttributeError:
+        print('Make sure to provide a Species object, w@nker')
+        raise