cannot reproduce end-to-end impact calc tutorial output

[dariomasante]

Describe the bug
Tutorial "end-to-end impact calc" does not output the shown value in jupyter block 13-14

To Reproduce
Steps to reproduce the behavior/error:
Run the tutorial up to block 13: https://github.com/CLIMADA-project/climada_python/blob/main/doc/tutorial/climada_engine_Impact.ipynb

Expected behavior
Example on on notebook reports 563366225 $, but I get 255124867 $

Additional context
Trying to do some comparisons, but cannot if it's not possible to reproduce the output of tutorials

[peanutfun]

@dariomasante Thank you for reporting this issue. I am assuming you refer to the average annual impact computed in block 13. The results present in the notebook are computed in January 2023, with the then latest version of Climada. There are several reasons why this number might have changed since then. The primary culprit might be the Hazard.frequency associated with the single events. Other than that, there might have been changes to the perturbation of cyclone tracks, the parametric wind field model, the exposure model, and the underlying data.

If the tests are passing on your system, I would not assume this is a bug, but rather that one or several parts of the modeling chain were updated, producing a different quantitative result.

[peanutfun]

More issues found in the review of the tutorial #1095

[luseverin]

I went back to climada v3.2.0 (f8b5c07), computed hazard, exposure and impact functions as well as impacts and I found that the discrepancies come from the hazards. It seems that the intensities of the tropical cyclone wind fields have changed between the current and the former version (somewhere before v3.2.0) that has been used to compute the tutorial. I have not yet been able to track the exact reason why, but I have found that 12 events over the 120 events have different intensities. Those are identified with the following ids:
array([ 30, 31, 34, 35, 38, 51, 58, 59, 98, 99, 116, 117])
I can also share some comparisons I have made between the current (top) vs old (bottom) tc events intensities:

In some cases the events seem markedly different (event_id=31), whereas in some other they seem rather similar (event_id=35). In some cases (event_id=34) it is not clear how they differ.
Anyone with a good overview on the recent changes in how we compute tropical cyclones' windfields has an idea of what might have happened?

[luseverin]

Possible cause for the differences in the windfields might come from updates to the ibtrack data.
TODO:

• Check if the differences indeed come from the changes in ibtrack tracks.

  climada_python/climada/hazard/tc_tracks.py

  Line 89 in bb95959

  IBTRACS_URL = (

[luseverin]

Switching the IBTracs data version between v04r01 (current: https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r01/access/netcdf) and v04r00 (climada 3.2.0: https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/netcdf/ ) does not seem to explain the large difference between the two versions as I obtain:

• climada 3.2.0 aai_agg = 561850498.0 $ with v04r01 vs aai_agg = 563366428.0 $ with v04r00
• climada 6.0.2.dev0 aai_agg = 257168070.0 $ with v04r01 vs aai_agg = 261368081.0 $with v04r00

So in my opinion this comes from changes in how we compute the windfields. What do you think?

[chahank]

Can you please check what is in TCTracks object for both versions? Because we changed some of the defaults in the get_ibtracks method.

[luseverin]

TCTracks objects look the same to me. The data has the same size and same fields. Did you have anything specific in mind?
Old data:

print(old_data.data[0])
<xarray.Dataset>
 Dimensions:                 (time: 349)
 Coordinates:
   * time                    (time) datetime64[ns] 2010-06-24T18:00:00 ... 201...
     lon                     (time) float64 -82.0 -82.01 -82.01 ... -101.8 -101.9
     lat                     (time) float64 15.9 15.91 15.92 ... 23.16 23.18 23.2
 Data variables:
     time_step               (time) float64 0.5 0.5 0.5 0.5 ... 0.5 0.5 0.5 0.5
     radius_max_wind         (time) float64 75.0 75.0 75.0 ... 47.33 46.17 45.0
     radius_oci              (time) float64 180.0 180.0 180.0 ... 130.0 120.0
     max_sustained_wind      (time) float64 25.0 25.0 25.0 ... 31.67 30.83 30.0
     central_pressure        (time) float64 1.007e+03 1.007e+03 ... 996.2 997.0
     environmental_pressure  (time) float64 1.011e+03 1.011e+03 ... 1.004e+03
     basin                   (time) <U2 'NA' 'NA' 'NA' 'NA' ... 'NA' 'NA' 'NA'
     on_land                 (time) bool False False False ... True True True
     dist_since_lf           (time) float64 nan nan nan nan ... 441.8 452.2 462.6
 Attributes:
     max_sustained_wind_unit:  kn
     central_pressure_unit:    mb
     name:                     ALEX
     sid:                      2010176N16278
     orig_event_flag:          1
     data_provider:            ibtracs_usa
     id_no:                    2010176016278.0
     category:                 2,

New data:

print(new_data.data[0])
<xarray.Dataset> Size: 31kB
 Dimensions:                 (time: 349)
 Coordinates:
   * time                    (time) datetime64[ns] 3kB 2010-06-24T18:00:00 ......
     lon                     (time) float64 3kB -82.0 -81.99 ... -101.8 -101.9
     lat                     (time) float64 3kB 15.9 15.89 15.88 ... 23.17 23.2
 Data variables:
     radius_max_wind         (time) float64 3kB 75.0 75.0 75.0 ... 46.33 45.0
     radius_oci              (time) float64 3kB 180.0 180.0 180.0 ... 130.0 120.0
     max_sustained_wind      (time) float64 3kB 25.0 25.0 25.0 ... 30.83 30.0
     central_pressure        (time) float64 3kB 1.007e+03 1.007e+03 ... 997.0
     environmental_pressure  (time) float64 3kB 1.011e+03 1.011e+03 ... 1.004e+03
     time_step               (time) float64 3kB 0.5 0.5 0.5 0.5 ... 0.5 0.5 0.5
     basin                   (time) <U2 3kB 'NA' 'NA' 'NA' ... 'NA' 'NA' 'NA'
     on_land                 (time) bool 349B False False False ... True True
     dist_since_lf           (time) float64 3kB nan nan nan ... 443.7 454.2 464.8
 Attributes:
     max_sustained_wind_unit:  kn
     central_pressure_unit:    mb
     orig_event_flag:          True
     data_provider:            ibtracs_usa
     category:                 2
     name:                     ALEX
     sid:                      2010176N16278
     id_no:                    2010176016278.0,

[chahank]

So, the track has the exact same coordinates and data variables? I wanted just to be sure that there is no difference in the underlying track data. Because as you diagnosed, we have an issue with the hazard, not the impact. So, we should look for the issue in the hazard. And here we have a hazard that is built from tracks. So, we should look at the differences in the data of the tracks and the resulting hazard objects. Not with plots though as this does not give that much information on what is really going on. Can one see which data are different in the two hazards? or in the tracks?

[peanutfun]

@luseverin I think you have to compare the dataset for each track, not just the first one. Only some of them appear to be different, as you already analyzed. xarray.testing.assert_equal might help here, as it reports well which values are different (if any)

[luseverin]

They do not have the same values of coordinates and variables, but I think this is expected as the underlying data from IBTrack is not the same between the two versions. When I use the two different TCTrack objects to generate the impacts, using CLIMADA 6.0.0, I obtain similar values for aai_agg (257168070.0 vs 259036421.0). So I still think it comes from changes in the windfield model.

[chahank]

So what we need to understand is why the updates to the tropical cyclone windfield model and/or to the tc from ibstracks lead to such large changes in the windfields. Hence, the impact should not be looked at anymore. Rather, we should compare the intensity matrices of the storms, and the xarray entries of the tracks.