new approach for plant DK

sample/dk_adv.py

@@ -64,6 +64,7 @@ Gebruik daarom de btA categorieen'''
 areas = hlp.gen_squares(x_ul=0, y_ul=625000, ncol=28, nrow=33, size=10000)
 bt_stats_prelim1 = rast.zonal_stats(vectors=areas, raster=bt_kaart.read(1), categorical=True, category_map=pxl2btA,
                                     affine=bt_kaart.affine)
+
 # Reclass btA categorieen naar btB categoriien
 bt_stats_prelim2 = []
 for stat in bt_stats_prelim1:

sample/dk_plant.py

@@ -10,19 +10,75 @@ Relation between species presence and vegetation types is sourced from DVN data,
 Hans Roelofsen, 23 June 2020
 """
 
-import os
+
+import itertools
 import pandas as pd
+from sample.mnp import helpers as hlp
+
 
 """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_mnp_dk = r'C:\apps\temp_geodata\MNP_DK\03_MNP_versie6_par_density_factors_BT2019_v2.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=',')
+mnp = pd.merge(left=mnp_sp_nms, right=mnp_sp_tax, how='left', left_on='Species_code', right_on='Species_code')
+mnp_dk = pd.read_csv(src_mnp_dk, sep=',', usecols=['Species_code', 'Land_type_code', 'Land_type_quality'])
 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
+'''Reduce sub-associations 00XX11a to 00XX11 and repair Beheertype codes in Input sheet vertaling_IPO_ASSO_2009.txt'''
+asso_x_bts.loc[:, 'syntaxon_gen'] = asso_x_bts.short.str.slice(stop=6).str.upper()
+asso_x_bts.loc[:, 'BT'] = asso_x_bts.Bc.apply(hlp.fix_bt)
+
+'''Reduce sub-associations in Input sheet q_syntaxa_soorten_frequentie_trouwsoorten.txt'''
+asso_x_sps.loc[:, 'syntaxon_gen'] = asso_x_sps.Syntaxon_code.str.slice(stop=6).str.upper()
+
+'''Reduce neergeschaalde beheertypen in MNP reference draagkrachten tabel'''
+sel = mnp_dk.loc[mnp_dk.Land_type_code.str.slice(start=7) == '00', :].index
+mnp_dk.loc[:, 'BT'] = None
+mnp_dk.loc[sel, 'BT'] = mnp_dk.loc[sel, 'Land_type_code'].str.slice(stop=6)
+
+'''Mappings between i) Associaties and Beheertypen (n:n) ii) name and species code'''
+bt2asso = {}
+for bt in list(set(asso_x_bts.BT)):
+    bt2asso[bt] = list(set(asso_x_bts.loc[asso_x_bts.BT == bt, 'syntaxon_gen']))
+name2code = dict(zip(mnp_sp_nms.Scientific_name, mnp_sp_nms.Species_code))
+
+'''New dataframe with all species <-> Beheertype combinations.'''
+dk = pd.DataFrame(data=[(a, b) for (a, b) in itertools.product(mnp.loc[mnp.Taxon_group == 'P', 'Scientific_name'],
+                                                               set(asso_x_bts.BT))], columns=['sp_name', 'BT'])
+dk.loc[:, 'sp_code'] = dk.sp_name.map(name2code)
+
+'''Summarize Frequentie for each BT <> Species combination'''
+qualitis = []
+mnp_ref = []
+for row in dk.itertuples():
+    # landtype_quality = gemiddelde frequentie van de soort in de associaties behorende bij het gevraagde BT
+    sel = asso_x_sps.loc[(asso_x_sps.Soortnaam == row.sp_name) & (asso_x_sps.syntaxon_gen.isin(bt2asso[row.BT])), :]
+    mean_freq = sel.Frequentie.mean()
+    qualitis.append(mean_freq)
+
+    # dit kan vast simpeler met groupby of zo, fuck it
+
+    if sel.shape[0] > 0:
+        print('Found {0} records for species {1} and BT {2}, with mean '
+              'Frequentie: {3}'.format(sel.shape[0], row.sp_name, row.BT, mean_freq))
+    mnp_ref.append(mnp_dk.loc[(mnp_dk.Species_code == row.sp_code) & (mnp_dk.BT == row.BT), 'Land_type_quality'].mean())
+
+dk.loc[:, 'land_type_q'] = qualitis
+dk.loc[:, 'mnp_land_type_q'] = mnp_ref
+dk.loc[dk.land_type_q.notna(), :].to_clipboard(sep=';', index=False)
+
+
+'''5. Reporting'''
+sp_done = len(set(dk.loc[dk.land_type_q.notna(), 'sp_name']))
+bt_done = len(set(dk.loc[dk.land_type_q.notna(), 'BT']))
+sp_all = len(set(dk.sp_name))
+u_bt = len(set(asso_x_bts.BT))  # nr of unique Beheertypen in WW's excel
+u_sp = len(set(asso_x_sps.Soortnaam))  # nr of unique species in WW's excel
+r_sp = len(set(mnp_sp_nms.Scientific_name).intersection(set(asso_x_sps.Soortnaam)))

sample/mnp/helpers.py

@@ -8,6 +8,7 @@ import geopandas as gp
 import pandas as pd
 import shapely
 import affine
+import numbers
 from datetime import datetime as dt
 from matplotlib import pyplot as plt
 
@@ -182,6 +183,19 @@ def gen_squares(x_ul, y_ul, nrow, ncol, size):
     return gp.GeoDataFrame(crs={"init": 'epsg:28992'}, geometry=shapes, data={'ID': ids, 'size': str(size)})
 
 
+def fix_bt(code_in):
+    '''Repair abbreviated BT codes'''
+    if isinstance(code_in, numbers.Number):
+        code_in = str(code_in)
+
+    head, tail = code_in.split('.')
+    if len(head) == 1:
+        pre = 'N0'
+    elif len(head) == 2:
+        pre = 'N'
+    else:
+        raise Exception('Unexpected BT code: {}'.format(code_in))
+    return '{0}{1}.{2}'.format(pre, head, tail)