Added feature to write_sera_to_bdg to narrow down com building classes

GDE_TOOLS_general.py

 """
-Copyright (C) 2020
-Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
+Copyright (C) 2020-2021
+  Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU Affero General Public License as
@@ -320,10 +320,10 @@ def write_sera_to_bdg(in_osm_id, open_fle_part, sera_dict, adm_level_from_OBM):
     Otherwise, he attribute "'Sto_SERA_inconsist_'+disag_case" is set to False. TO DO: IN THE FUTURE WE MIGHT
     WANT TO DO SOMETHING DIFFERENT WITH CASES LIKE THIS. NEED TO DISCUSS.
     """
-    err_str= ''
+    err_str = ''
     for disag_case in sera_dict.keys():
-        if len(sera_dict[disag_case].keys())!=0: # SERA HDF5 contains info for this cell and admin unit         
-            if sera_dict[disag_case]['Adm_Level']!=adm_level_from_OBM:
+        if len(sera_dict[disag_case].keys()) != 0:  # SERA HDF5 contains info for this cell and admin unit         
+            if sera_dict[disag_case]['Adm_Level'] != adm_level_from_OBM:
                 err_str= 'ERROR IN ADM LEVEL OF SERA HDF5 FILES VS OBM DATABASE FOR '+disag_case+'.'
             else:
                 if 'SERA_classes'+'_'+disag_case in list(open_fle_part[str(in_osm_id)]):
@@ -331,14 +331,33 @@ def write_sera_to_bdg(in_osm_id, open_fle_part, sera_dict, adm_level_from_OBM):
                 if 'SERA_vals'+'_'+disag_case in list(open_fle_part[str(in_osm_id)]):
                     del open_fle_part[str(in_osm_id)]['SERA_vals'+'_'+disag_case]
                 # Narrowing down the distribution based on the number of storeys, only if available:
-                num_sto= open_fle_part[str(in_osm_id)].attrs['Storeys']
-                keep_class= np.array([True for i in range(0, len(sera_dict[disag_case]['SERA_classes']))]) # keep all classes unles removed
-                if num_sto!='UNK':
-                    if float(num_sto)>0.1:
-                        num_sto_lower, num_sto_upper= gdet_sera.interpret_num_storeys_from_building_class(sera_dict[disag_case]['SERA_classes'])
+                num_sto = open_fle_part[str(in_osm_id)].attrs['Storeys']
+                keep_class_sto = np.array([True for i in range(0, len(sera_dict[disag_case]['SERA_classes']))])  # keep all classes unless removed
+                if num_sto != 'UNK':
+                    if float(num_sto) > 0.1:
+                        num_sto_lower, num_sto_upper = gdet_sera.interpret_num_storeys_from_building_class(sera_dict[disag_case]['SERA_classes'])
                         for k, theo_class in enumerate(sera_dict[disag_case]['SERA_classes']):                   
-                            if not np.logical_and(np.ceil(float(num_sto))>=num_sto_lower[k], np.ceil(float(num_sto))<=num_sto_upper[k]): # np.ceil(float), not int, because there are OSM buildings with non-integer number of storeys (ceiling needed bec if e.g. num_sto='4.5' and the bdg classes I am comparing against have H:4 and H:5, the code would not assign it to any of the two classes otherwise)
-                                keep_class[k]= False # remove class  
+                            if not np.logical_and(np.ceil(float(num_sto)) >= num_sto_lower[k], np.ceil(float(num_sto)) <= num_sto_upper[k]):  # np.ceil(float), not int, because there are OSM buildings with non-integer number of storeys (ceiling needed bec if e.g. num_sto='4.5' and the bdg classes I am comparing against have H:4 and H:5, the code would not assign it to any of the two classes otherwise)
+                                keep_class_sto[k] = False # remove class  
+                # Narrowing down the distribution based on commercial occupancies:
+                occupancy = open_fle_part[str(in_osm_id)].attrs["Occup"]
+                keep_class_com = np.array([True for i in range(0, len(sera_dict[disag_case]['SERA_classes']))])  # keep all classes unless removed
+                if (occupancy == "RES3"  # Hotels
+                    or occupancy == "COM5"):  # Restaurants
+                    for k, theo_class in enumerate(sera_dict[disag_case]['SERA_classes']):
+                        if "Hotels" not in theo_class:  # In SERA, "Hotels" includes hotels and restaurants
+                            keep_class_com[k] = False # remove class
+                elif occupancy == "COM3":  # Offices
+                    for k, theo_class in enumerate(sera_dict[disag_case]['SERA_classes']):
+                        if "Offices" not in theo_class:
+                            keep_class_com[k] = False # remove class
+                elif (occupancy == "COM1"  # Retail trade
+                      or occupancy == "COM2"):  # Wholesale trade and storage (warehouse)
+                    for k, theo_class in enumerate(sera_dict[disag_case]['SERA_classes']):
+                        if "Trade" not in theo_class:  # In SERA, "Trade" includes wholesale and retail trade
+                            keep_class_com[k] = False # remove class
+                # Combine results from both types of narrowing down cases:
+                keep_class = np.logical_and(keep_class_sto, keep_class_com)
                 if keep_class.sum()>0.1: # i.e. if there is at least one class left as "to keep", i.e. at least one "True" in keep_class
                     # Write the classes and values to keep (it is all of them if num_sto is unknown):
                     dt=h5py.special_dtype(vlen=str) # examples on using special_dtype: https://www.programcreek.com/python/example/105243/h5py.special_dtype