I'm getting negative 50-year impacts for tropical storm winds when I calculate them at centroid level.

I'm using the Impact.local_exceedance_imp method to calculate impacts by return period and centroid. The method takes centroids and your desired return periods, and for each centroid estimates the retunr period impacts with the Impact._cen_return_imp method.

My understanding of Impact._cen_return_imp is that it takes the ranked impacts and cumulative event frequencies for the centroid and calculates the line of best fit through them (taking log of the frequencies first). It then estimates the impact for the desired return periods from this line (I've copied the code for the method below).

This method is different to the Impact.def_calc_freq method (used to calculate impacts by return periods for all centroids). In contrast the estimations here use np.interp to interpolate between the calculated fequencies and impacts.

I think I understand the reasoning for the difference: a single centroid might not have many events, so fitting a logarithmic curve to the data instead of interpolating might avoid some weird results?

I want to discuss it for two reasons

1. It's unexpected. I would expect to set a different method of calculation (especially modelled through best fit) explicitly.
2. Negative impacts for lower return periods in my analysis (even 50 years) is bad.

Thoughts? Does anyone know more about the history of the method?

EDIT: same question applies to Hazard and Hazard._cen_return_inten

    @staticmethod
    def _cen_return_imp(imp, freq, imp_th, return_periods):
        """From ordered impact and cummulative frequency at centroid, get
        exceedance impact at input return periods.
        Parameters:
            imp (np.array): sorted impact at centroid
            freq (np.array): cummulative frequency at centroid
            imp_th (float): impact threshold
            return_periods (np.array): return periods
        Returns:
            np.array
        """
        imp_th = np.asarray(imp > imp_th).squeeze()
        imp_cen = imp[imp_th]
        freq_cen = freq[imp_th]
        if not imp_cen.size:
            return np.zeros((return_periods.size,))
        try:
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                pol_coef = np.polyfit(np.log(freq_cen), imp_cen, deg=1)
        except ValueError:
            pol_coef = np.polyfit(np.log(freq_cen), imp_cen, deg=0)
        imp_fit = np.polyval(pol_coef, np.log(1 / return_periods))
        wrong_inten = (return_periods > np.max(1 / freq_cen)) & np.isnan(imp_fit)
        imp_fit[wrong_inten] = 0.

        return imp_fit