Feature/reassign centr bydflt

climada/engine/calibration_opt.py

@@ -66,26 +66,26 @@ def calib_instance(hazard, exposure, impact_func, df_out=pd.DataFrame(),
             DataFrame with modelled impact written to rows for each year
             or event.
     """
-    IFS = ImpactFuncSet()
-    IFS.append(impact_func)
+    ifs = ImpactFuncSet()
+    ifs.append(impact_func)
     impacts = Impact()
-    impacts.calc(exposure, IFS, hazard)
+    impacts.calc(exposure, ifs, hazard, assign_centroids=False)
     if yearly_impact:  # impact per year
-        IYS = impacts.calc_impact_year_set(all_years=True)
+        iys = impacts.calc_impact_year_set(all_years=True)
         # Loop over whole year range:
         if df_out.empty | df_out.index.shape[0] == 1:
-            for cnt_, year in enumerate(np.sort(list((IYS.keys())))):
+            for cnt_, year in enumerate(np.sort(list((iys.keys())))):
                 if cnt_ > 0:
                     df_out.loc[cnt_] = df_out.loc[0]  # copy info from first row
-                if year in IYS:
-                    df_out.loc[cnt_, 'impact_CLIMADA'] = IYS[year]
+                if year in iys:
+                    df_out.loc[cnt_, 'impact_CLIMADA'] = iys[year]
                 else:
                     df_out.loc[cnt_, 'impact_CLIMADA'] = 0.0
                 df_out.loc[cnt_, 'year'] = year
         else:
-            years_in_common = df_out.loc[df_out['year'].isin(np.sort(list((IYS.keys())))), 'year']
+            years_in_common = df_out.loc[df_out['year'].isin(np.sort(list((iys.keys())))), 'year']
             for cnt_, year in years_in_common.iteritems():
-                df_out.loc[df_out['year'] == year, 'impact_CLIMADA'] = IYS[year]
+                df_out.loc[df_out['year'] == year, 'impact_CLIMADA'] = iys[year]
 
 
     else:  # impact per event
@@ -138,21 +138,21 @@ def init_impf(impf_name_or_instance, param_dict, df_out=pd.DataFrame(index=[0]))
         (index=0) defined with values. The impact function parameters from
         param_dict are represented here.
     """
-    ImpactFunc_final = None
+    impact_func_final = None
     if isinstance(impf_name_or_instance, str):
         if impf_name_or_instance == 'emanuel':
-            ImpactFunc_final = ImpfTropCyclone.from_emanuel_usa(**param_dict)
-            ImpactFunc_final.haz_type = 'TC'
-            ImpactFunc_final.id = 1
+            impact_func_final = ImpfTropCyclone.from_emanuel_usa(**param_dict)
+            impact_func_final.haz_type = 'TC'
+            impact_func_final.id = 1
             df_out['impact_function'] = impf_name_or_instance
     elif isinstance(impf_name_or_instance, impact_funcs.ImpactFunc):
-        ImpactFunc_final = change_impf(impf_name_or_instance, param_dict)
+        impact_func_final = change_impf(impf_name_or_instance, param_dict)
         df_out['impact_function'] = ('given_' +
-                                     ImpactFunc_final.haz_type +
-                                     str(ImpactFunc_final.id))
+                                     impact_func_final.haz_type +
+                                     str(impact_func_final.id))
     for key, val in param_dict.items():
         df_out[key] = val
-    return ImpactFunc_final, df_out
+    return impact_func_final, df_out
 
 def change_impf(impf_instance, param_dict):
     """apply a shifting or a scaling defined in param_dict to the impact
@@ -333,6 +333,7 @@ def calib_all(hazard, exposure, impf_name_or_instance, param_full_dict,
     # prepare hazard and exposure
     region_ids = list(np.unique(exposure.region_id))
     hazard_type = hazard.tag.haz_type
+    exposure.assign_centroids(hazard)
     # prepare impact data
     if isinstance(impact_data_source, pd.DataFrame):
         df_impact_data = impact_data_source
@@ -347,8 +348,8 @@ def calib_all(hazard, exposure, impf_name_or_instance, param_full_dict,
     for param_dict in params_generator:
         print(param_dict)
         df_out = copy.deepcopy(df_impact_data)
-        ImpactFunc_final, df_out = init_impf(impf_name_or_instance, param_dict, df_out)
-        df_out = calib_instance(hazard, exposure, ImpactFunc_final, df_out, yearly_impact)
+        impact_func_final, df_out = init_impf(impf_name_or_instance, param_dict, df_out)
+        df_out = calib_instance(hazard, exposure, impact_func_final, df_out, yearly_impact)
         if df_result is None:
             df_result = copy.deepcopy(df_out)
         else:
@@ -398,6 +399,7 @@ def calib_optimize(hazard, exposure, impf_name_or_instance, param_dict,
     # prepare hazard and exposure
     region_ids = list(np.unique(exposure.region_id))
     hazard_type = hazard.tag.haz_type
+    exposure.assign_centroids(hazard)
     # prepare impact data
     if isinstance(impact_data_source, pd.DataFrame):
         df_impact_data = impact_data_source
@@ -408,8 +410,8 @@ def calib_optimize(hazard, exposure, impf_name_or_instance, param_dict,
         else:
             raise ValueError('other impact data sources not yet implemented.')
     # definie specific function to
-    def specific_calib(x):
-        param_dict_temp = dict(zip(param_dict.keys(), x))
+    def specific_calib(values):
+        param_dict_temp = dict(zip(param_dict.keys(), values))
         print(param_dict_temp)
         return calib_instance(hazard, exposure,
                               init_impf(impf_name_or_instance, param_dict_temp)[0],
@@ -427,8 +429,8 @@ def specific_calib(x):
                 {'type': 'ineq', 'fun': lambda x: x[1]}]
 
 
-    x0 = list(param_dict.values())
-    res = minimize(specific_calib, x0,
+    values = list(param_dict.values())
+    res = minimize(specific_calib, values,
                    # bounds=bounds,
                    # bounds=((0.0, np.inf), (0.0, np.inf), (0.0, 1.0)),
                    constraints=cons,

climada/engine/cost_benefit.py

@@ -158,7 +158,7 @@ def __init__(self):
         self.imp_meas_present = dict()
 
     def calc(self, hazard, entity, haz_future=None, ent_future=None,
-             future_year=None, risk_func=risk_aai_agg, imp_time_depen=None, save_imp=False):
+             future_year=None, risk_func=risk_aai_agg, imp_time_depen=None, save_imp=False, assign_centroids=True):
         """Compute cost-benefit ratio for every measure provided current
         and, optionally, future conditions. Present and future measures need
         to have the same name. The measures costs need to be discounted by the user.
@@ -203,6 +203,14 @@ def calc(self, hazard, entity, haz_future=None, ent_future=None,
         if future_year is None and ent_future is None:
             future_year = entity.exposures.ref_year
 
+        # assign centroids
+        if assign_centroids:
+            entity.exposures.assign_centroids(hazard, overwrite=True)
+            if ent_future:
+                ent_future.exposures.assign_centroids(
+                    haz_future if haz_future else hazard, overwrite=True
+                )
+
         if not haz_future and not ent_future:
             self.future_year = future_year
             self._calc_impact_measures(hazard, entity.exposures,
@@ -525,11 +533,11 @@ def plot_waterfall(hazard, entity, haz_future, ent_future,
         future_year = ent_future.exposures.ref_year
 
         imp = Impact()
-        imp.calc(entity.exposures, entity.impact_funcs, hazard)
+        imp.calc(entity.exposures, entity.impact_funcs, hazard, assign_centroids=False)
         curr_risk = risk_func(imp)
 
         imp = Impact()
-        imp.calc(ent_future.exposures, ent_future.impact_funcs, haz_future)
+        imp.calc(ent_future.exposures, ent_future.impact_funcs, haz_future, assign_centroids=False)
         fut_risk = risk_func(imp)
 
         if not axis:
@@ -542,7 +550,7 @@ def plot_waterfall(hazard, entity, haz_future, ent_future,
         # changing future
         # socio-economic dev
         imp = Impact()
-        imp.calc(ent_future.exposures, ent_future.impact_funcs, hazard)
+        imp.calc(ent_future.exposures, ent_future.impact_funcs, hazard, assign_centroids=False)
         risk_dev = risk_func(imp)
         LOGGER.info('Risk with development at {:d}: {:.3e}'.format(future_year, risk_dev))
 
@@ -687,7 +695,7 @@ def plot_waterfall_accumulated(self, hazard, entity, ent_future,
         time_dep = self._time_dependency_array(imp_time_depen)
         # socio-economic dev
         imp = Impact()
-        imp.calc(ent_future.exposures, ent_future.impact_funcs, hazard)
+        imp.calc(ent_future.exposures, ent_future.impact_funcs, hazard, assign_centroids=False)
         risk_dev = self._npv_unaverted_impact(risk_func(imp), entity.disc_rates,
                                               time_dep, curr_risk)
         LOGGER.info('Total risk with development at {:d}: {:.3e}'.format(
@@ -760,7 +768,7 @@ def _calc_impact_measures(self, hazard, exposures, meas_set, imp_fun_set,
         # compute impact without measures
         LOGGER.debug('%s impact with no measure.', when)
         imp_tmp = Impact()
-        imp_tmp.calc(exposures, imp_fun_set, hazard)
+        imp_tmp.calc(exposures, imp_fun_set, hazard, assign_centroids=False)
         impact_meas[NO_MEASURE] = dict()
         impact_meas[NO_MEASURE]['cost'] = (0, 0)
         impact_meas[NO_MEASURE]['risk'] = risk_func(imp_tmp)
@@ -772,7 +780,7 @@ def _calc_impact_measures(self, hazard, exposures, meas_set, imp_fun_set,
         # compute impact for each measure
         for measure in meas_set.get_measure(hazard.tag.haz_type):
             LOGGER.debug('%s impact of measure %s.', when, measure.name)
-            imp_tmp, risk_transf = measure.calc_impact(exposures, imp_fun_set, hazard)
+            imp_tmp, risk_transf = measure.calc_impact(exposures, imp_fun_set, hazard, assign_centroids=False)
             impact_meas[measure.name] = dict()
             impact_meas[measure.name]['cost'] = (measure.cost, measure.risk_transf_cost_factor)
             impact_meas[measure.name]['risk'] = risk_func(imp_tmp)

climada/engine/forecast.py

@@ -292,12 +292,12 @@ def calc(self, force_reassign=False):
             default is false.
         """
         # calc impact
+        if self.hazard:
+            self.exposure.assign_centroids(self.hazard[0], overwrite=force_reassign)
         for ind_i, haz_i in enumerate(self.hazard):
-            # force reassign
-            if force_reassign:
-                self.exposure.assign_centroids(haz_i)
             self._impact[ind_i].calc(
-                self.exposure, self.vulnerability, haz_i, save_mat=True
+                self.exposure, self.vulnerability, haz_i,
+                save_mat=True, assign_centroids=False
             )
 
     def plot_imp_map(

climada/engine/impact.py

@@ -184,7 +184,7 @@ def __init__(self,
 
 
 
-    def calc(self, exposures, impact_funcs, hazard, save_mat=False):
+    def calc(self, exposures, impact_funcs, hazard, save_mat=False, assign_centroids=True):
         """This function is deprecated, use ImpactCalc.impact
         and ImpactCalc.insured_impact instead.
         """
@@ -198,11 +198,11 @@ def calc(self, exposures, impact_funcs, hazard, save_mat=False):
                 " for insured impacts instead. For non-insured impacts "
                 "please use ImpactCalc().impact()"
                 )
-            self.__dict__ = impcalc.insured_impact(save_mat).__dict__
+            self.__dict__ = impcalc.insured_impact(save_mat, assign_centroids).__dict__
         else:
             LOGGER.warning("The use of Impact().calc() is deprecated. "
                            "Use ImpactCalc().impact() instead.")
-            self.__dict__ = impcalc.impact(save_mat).__dict__
+            self.__dict__ = impcalc.impact(save_mat, assign_centroids).__dict__
 
 #TODO: new name
     @classmethod
@@ -1027,9 +1027,11 @@ def video_direct_impact(exp, impf_set, haz_list, file_name='',
         imp_list = []
         exp_list = []
         imp_arr = np.zeros(len(exp.gdf))
+        # assign centroids once for all
+        exp.assign_centroids(haz_list[0])
         for i_time, _ in enumerate(haz_list):
             imp_tmp = Impact()
-            imp_tmp.calc(exp, impf_set, haz_list[i_time])
+            imp_tmp.calc(exp, impf_set, haz_list[i_time], assign_centroids=False)
             imp_arr = np.maximum(imp_arr, imp_tmp.eai_exp)
             # remove not impacted exposures
             save_exp = imp_arr > imp_thresh

climada/engine/impact_calc.py

@@ -61,7 +61,8 @@ def __init__(self,
         self.exposures = exposures
         self.impfset = impfset
         self.hazard = hazard
-        self._orig_exp_idx = np.arange(self.exposures.gdf.shape[0]) #exposures index to use for matrix reconstruction
+        # exposures index to use for matrix reconstruction
+        self._orig_exp_idx = np.arange(self.exposures.gdf.shape[0])
 
     @property
     def n_exp_pnt(self):
@@ -101,13 +102,21 @@ def cover(self):
         if 'cover' in self.exposures.gdf.columns:
             return self.exposures.gdf['cover'].to_numpy()
 
-    def impact(self, save_mat=True):
+    def impact(self, save_mat=True, assign_centroids=True):
         """Compute the impact of a hazard on exposures.
 
         Parameters
         ----------
-        save_mat : bool
+        save_mat : bool, optional
             if true, save the total impact matrix (events x exposures)
+            Default: True
+        assign_centroids : bool, optional
+            indicates whether centroids are re-assigned to the self.exposures object
+            or kept from previous impact calculation with a hazard of the same hazard type.
+            Centroids assignment is an expensive operation; set this to ``False`` to save
+            computation time if the centroids have not changed since the last impact
+            calculation.
+            Default: True
 
         Examples
         --------
@@ -126,7 +135,7 @@ def impact(self, save_mat=True):
         the column is added to the exposures geodataframe.
         """
         impf_col = self.exposures.get_impf_column(self.hazard.haz_type)
-        exp_gdf = self.minimal_exp_gdf(impf_col)
+        exp_gdf = self.minimal_exp_gdf(impf_col, assign_centroids)
         if exp_gdf.size == 0:
             return self._return_empty(save_mat)
         LOGGER.info('Calculating impact for %s assets (>0) and %s events.',
@@ -136,7 +145,7 @@ def impact(self, save_mat=True):
 
 #TODO: make a better impact matrix generator for insured impacts when
 # the impact matrix is already present
-    def insured_impact(self, save_mat=False):
+    def insured_impact(self, save_mat=False, assign_centroids=True):
         """Compute the impact of a hazard on exposures with a deductible and/or
         cover.
 
@@ -147,6 +156,13 @@ def insured_impact(self, save_mat=False):
         ----------
         save_mat : bool
             if true, save the total impact matrix (events x exposures)
+        assign_centroids : bool, optional
+            indicates whether centroids are re-assigned to the self.exposures object
+            or kept from previous impact calculation with a hazard of the same hazard type.
+            Centroids assignment is an expensive operation; set this to ``False`` to save
+            computation time if the centroids have not changed since the last impact
+            calculation.
+            Default: True
 
         Examples
         --------
@@ -169,7 +185,7 @@ def insured_impact(self, save_mat=False):
                                  "Please set exposures.gdf.cover"
                                  "and/or exposures.gdf.deductible")
         impf_col = self.exposures.get_impf_column(self.hazard.haz_type)
-        exp_gdf = self.minimal_exp_gdf(impf_col)
+        exp_gdf = self.minimal_exp_gdf(impf_col, assign_centroids)
         if exp_gdf.size == 0:
             return self._return_empty(save_mat)
         LOGGER.info('Calculating impact for %s assets (>0) and %s events.',
@@ -238,18 +254,29 @@ def _return_empty(self, save_mat):
         return Impact.from_eih(self.exposures, self.impfset, self.hazard,
                         at_event, eai_exp, aai_agg, imp_mat)
 
-    def minimal_exp_gdf(self, impf_col):
+    def minimal_exp_gdf(self, impf_col, assign_centroids):
         """Get minimal exposures geodataframe for impact computation
 
         Parameters
         ----------
         exposures : climada.entity.Exposures
         hazard : climada.Hazard
-        impf_col: str
-            name of the impact function column in exposures.gdf
-
-        """
-        self.exposures.assign_centroids(self.hazard, overwrite=False)
+        impf_col : str
+            Name of the impact function column in exposures.gdf
+        assign_centroids : bool
+            Indicates whether centroids are re-assigned to the self.exposures object
+            or kept from previous impact calculation with a hazard of the same hazard type.
+            Centroids assignment is an expensive operation; set this to ``False`` to save
+            computation time if the centroids have not changed since the last impact
+            calculation.
+        """
+        if assign_centroids:
+            self.exposures.assign_centroids(self.hazard, overwrite=True)
+        elif self.hazard.centr_exp_col not in self.exposures.gdf.columns:
+            raise ValueError("'assign_centroids' is set to 'False' but no centroids are assigned"
+                             f" for the given hazard type ({self.hazard.tag.haz_type})."
+                             " Run 'exposures.assign_centroids()' beforehand or set"
+                             " 'assign_centroids' to 'True'")
         mask = (
             (self.exposures.gdf.value.values != 0)
             & (self.exposures.gdf[self.hazard.centr_exp_col].values >= 0)
@@ -260,7 +287,8 @@ def minimal_exp_gdf(self, impf_col):
             )
         if exp_gdf.size == 0:
             LOGGER.warning("No exposures with value >0 in the vicinity of the hazard.")
-        self._orig_exp_idx = mask.nonzero()[0]  #update index of kept exposures points in exp_gdf within the full exposures
+        self._orig_exp_idx = mask.nonzero()[0]  # update index of kept exposures points in exp_gdf
+                                                # within the full exposures
         return exp_gdf
 
     def imp_mat_gen(self, exp_gdf, impf_col):
@@ -381,10 +409,11 @@ def impact_matrix(self, exp_values, cent_idx, impf):
         scipy.sparse.csr_matrix
             Impact per event (rows) per exposure point (columns)
         """
-        n_exp_pnt = len(cent_idx) #implicitly checks in matrix assignement whether len(cent_idx) == len(exp_values)
+        n_exp_pnt = len(cent_idx)  # implicitly checks in matrix assignement whether
+                                   # len(cent_idx) == len(exp_values)
         mdr = self.hazard.get_mdr(cent_idx, impf)
         fract = self.hazard.get_fraction(cent_idx)
-        exp_values_csr = sparse.csr_matrix( #vector 1 x exp_size
+        exp_values_csr = sparse.csr_matrix(  # vector 1 x exp_size
             (exp_values, np.arange(n_exp_pnt), [0, n_exp_pnt]),
             shape=(1, n_exp_pnt))
         return fract.multiply(mdr).multiply(exp_values_csr)
@@ -393,7 +422,9 @@ def stitch_impact_matrix(self, imp_mat_gen):
         """
         Make an impact matrix from an impact sub-matrix generator
         """
-        data, row, col = np.hstack([  #rows=events index, cols=exposure point index within self.exposures
+        # rows: events index
+        # cols: exposure point index within self.exposures
+        data, row, col = np.hstack([
             (mat.data, mat.nonzero()[0], self._orig_exp_idx[idx][mat.nonzero()[1]])
             for mat, idx in imp_mat_gen
             ])