Negative TC impacts: a question about how we calculate exceedences for centroids

[ChrisFairless]

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

[ChrisFairless]

I've reimplemented this with the same linear interpolation method as Impact.def_calc_freq on the branch feature/impact-alt-exceedence-method, just for my own use for now.

[chahank]

Great that you bring this up! I was just recently having similar troubles with this function.

The passage to the log-space is probably sensible due to the scarse and multi-scale data. However, the fit should be done in a more sophisticated way to avoid negative values. How about forcing the fit to pass in the first point?

[chahank]

To extend a bit: I agree that Hazard._cen_return_imp() and Impact.calc_freq_curve() work in a fundamentally different way.

Impact.calc_freq_curve(): only linear interpolation between all the defined points Hazard._cen_return_imp(): linear or constant fitting of the curve in log space

I think your suggestion to simply use the interpolation for Hazard._cen_return_imp() could be a quick and dirty fix. Could you do some testing to get a feeling for how good this solution is and whether a more sophisticated fitting would be needed?

[ChrisFairless]

It's a difficult question. I wonder if the assumption that exceedance curves are log-linear is too strong. A lot of the curves generated by calc_freq_curve don't look log-linear (this might be another topic of conversation).

My opinion is that we shouldn't be fitting a model that makes such strong assumptions without the user's knowledge, and that maybe the user shouldn't expect CLIMADA to make a good exceedance curve from insufficient data.

However if we want to stay with the log-linear interpolation, I think forcing the fit to pass through the first point (lowest RP, lowest impact) is a really good idea, since that point won't be very sensitive.

Another option is that we interpolate as in calc_freq_curve, but log-linear, and extrapolating beyond the highest RP in the data so that we get higher impacts for unrepresented return periods without anything ridiculous happening. This would change the results very little from calc_freq_curves linear interpolation, except at high return periods and when very few events are available. But again, this is not really a solution for an event set that doesn't cover low-frequency events well at a single location.

[ChrisFairless]

Let me know how we might test interpolation as a quick and dirty fix. I've only checked it in locations with good event coverage, where it works fine. (Maybe we subsample the event set and check for errors at different return periods against the full 'truth' event set?)

[davidnbresch]

exceedance probability curves follow a GEV (for lack of a better reference: https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution), usually a Pareto-curve is a good aproximation (and has only one free parameter, alpha, https://en.wikipedia.org/wiki/Pareto_distribution), sometimes the specila forms like Fréchet, Weibull and Gumbel versions are fine, too.

We shall always show the original, unsmoothed curve, but any further math/approximation/smoothing/fitting is welcome.

One can imagine impacts that would not follow a GEV and for sure ones where no smoothing or other operations is necessary.

On 19 Apr 2021, at 14:27, Chris Fairless @.***> wrote:

Let me know how we might test interpolation as a quick and dirty fix. I've only checked it in locations with good event coverage, where it works fine. (Maybe we subsample the event set and check for errors at different return periods against the full 'truth' event set?)

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/CLIMADA-project/climada_python/issues/209#issuecomment-822427024, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABW3G5HL46ZZLETQWTXB6OTTJQOS7ANCNFSM427N2VQA.

[chahank]

@davidnbresch : the problem with showing the the original, unsmoothed curve for single centroid points is that some points will have very few events. Thus, the curve is not really a curve and we have to decide how to show/compute it. @ChrisFairless proposed two solutions, and we could add a 3rd one:

• linear interpolation (such as is done for the curve for all centroids)
• log-linear interpolation (maybe better coverage of large events)
• Pareto-curve fitting (better overall function, however contains a model assumption)
• (Bad) currently, it is a log-linear curve fitting (leads to negative loss values among others)

[chahank]

Let me know how we might test interpolation as a quick and dirty fix. I've only checked it in locations with good event coverage, where it works fine. (Maybe we subsample the event set and check for errors at different return periods against the full 'truth' event set?)

I am not sure whether this will work right? The whole point is to see how the curves difffer locally from the global average no?

[ChrisFairless]

I am not sure whether this will work right? The whole point is to see how the curves difffer locally from the global average no?

I think we can't test against a (spatially) 'global' average due to the different spatial scales involved (if I've understood you correctly). The exceedence curve for a location will always look different to the exceedence calculated for a larger region around it, because of how the numbers are aggregated (e.g. the 50-year max wind for Miami is a different distribution than the 50-year max wind for Florida).

So my experiment might go like this:

• Choose some locations/regions where we have many events.
• For each of those locations, repeatedly sample n events, infer an exceedence curve using each of the methods above, and calculate an average error.
• See how the average error changes with method and with n

We should also repeat the experiment to see how well the methods perform when they're trying to approximate curves from an ideal Weibull distribution (and also check how far our chosen locations are from a Weibull distribution...)

[chahank]

That sounds very reasonable, although it will be difficult to make a global 'calibration' in this way. Probably one should repeat the experiment at least for all the major basins separately.

[micgloor]

@chahank Would you have any updates on this issue? Thanks a lot & best regards, Michael

[ChrisFairless]

Hey, we've just run into this problem again while working on European wind storms. Luca spent a week looking for bugs that were giving weird 10-year return period losses, and we eventually realised it was this issue.

We fixed it using the feature/impact-alt-exceedence-method branch.

But yeah we need to agree on a fix, this unexpected behaviour isn't good!

[chahank]

Well, having a local approximation of an exceedance curve is a shaky exercise. But, please, propose an update. I do not think this is a bug-fix, rather a method update.

[ChrisFairless]

Ok, here's my suggestion, building on the earlier points:

We add an extra parameter to exposure-level and centroid-level return period calculations called e.g. curve_fit and create a util file with different functions that can be to fit curves, including the current log-linear fit.

We set the default value to be a linear fit function, because that matches the behaviour elsewhere in CLIMADA.

Users can add their own methods as they need them (e.g. Pareto, log-linear interpolation).

If we do this I'm concerned about backward compatibility. This has the potential to change a lot of existing analyses, potentially by large amounts. I think they will be improved, but there won't be any indication to the user what's changed. We have options:

• Do nothing (easiest)
• Add a warning to inform people and get rid of it in a year. CLIMADA loves warning you about things
• Implement the changes as a new method in the new ImpactCalc class instead and don't change the methods in Impact (I don't like this one)

I'll start working on this unless I hear otherwise.

[chahank]

@ChrisFairless : thanks for the suggestion!

• What do you mean by centroid-level return period calculations? Do you mean for Hazard objects? Otherwise, it is always Exposure points as the impact is only defined at exposure points, and not at centroids.
• Probably it makes the most sense to use interpolation methods from an established package such as scipy.
• In general, the return period exceedance curve is not fitted, so the default is not really linear. Adding a fit on top of the non-fitted data is something that would not break back compatibility if done well. Basically, always compute the values at the discrete points available. Then, in addition, provide a fitting method for visualization and further analysis.
• On possibility could be to make more serious use of the object ImpactFreqCurve. This object could have a fit method, which returns a fit to the data based on different methods. These returned values can be either be stored, or used directly for further analysis and plotting.

[ChrisFairless]

In order:

• Centroid-level vs Exposure-level: yes – I meant for for Hazard return period calculations as well - we should use the same approach for both. Right now the centroid-level calculations are done by linear interpolation, in contrast to the exposure-level calculations.
• Yes, let's use scipy etc
• Ah I've been super unclear with the language here. When I've been saying 'linear' I mean piecewise linear interpolation between points from the raw hazard/impact data (i.e. calculated from cumulative frequency and modelled hazard/impact). That's in contrast to the curve that's causing the negative impacts which is log-linear regression.
• I like the idea of using the ImpactFreqCurve. I can look into this.

[micgloor]

Dear @chahank and @ChrisFairless, I was wondering if you agreed on a way forward to fix this issue or bring Chris' fix to the master branch?

Thanks a lot for the update!

Best regards Michael

[chahank]

Dear @micgloor ,

unfortunately no, this fix is not yet ready for merging, and @ChrisFairless is not working on it. Would you like to take this over? That would be very helpful. Also, would you also be so kind and let us know who you are?

[micgloor]

Dear @chahank

Sorry for not identifying myself ;-) I though you might have recognized it from the username. I'm Michael Gloor from Correntics and we met a while ago at ETH. We're currently collaborating with @davidnbresch and @aleeciu in the NCCS project.

As this fix will likely affect the reproducibility of previous results and I'm not sure if I'm in the best position to assess the implications on other projects, I think it would be good if the CLIMADA core development team could propose the way forward. Based on the previous discussion above, there seem to be different options.

Thanks & best regards Michael

[chahank]

Dear Michael, thanks for the answer and the context!

For the fix, there is currently no one working on this as far as I know. The original developer @ChrisFairless is as I understand not planning to do it in the near future. We would very much welcome it if you can propose an implementation based on the above discussion in the form of a pull request.