Tropical Cyclone Image Artifact Removal/Smoothing

[haicyu]

Is your feature request related to a problem? Please describe.
A clear and concise description of what the problem is. Ex. I'm always frustrated when [...]
Spurious image artifacts are making the hazard intensity profiles for tropical cyclones unusable at a high resolution.
Tropical cyclone storm tracks are printed into hazard intensity profiles at high spatial resolution with low temporal resolution. This results in spatial artifacts instead of a smooth profile of a tropical cyclone traveling through space over time. Is there any way to plot a smoother intensity profile for a given tropical cyclone?

Describe the solution you'd like
A clear and concise description of what you want to happen.

Describe alternatives you've considered
A clear and concise description of any alternative solutions or features you've considered.
I've increased the number of synthetic tracks generated to get more smoothness by having them overlap with each other, but I've reached the limits of my system. Currently my model takes 4 hours to build, and I am still seeing these image artifacts. Is there any other way to get some smoothing? We should get a smooth intensity profile for a tropical cyclone along the direction of tropical cyclone travel and in a general region.

To Reproduce
Steps to reproduce the behavior/error:

1. Run CLIMADA at high resolution, with spacing between centroids less than 0.1 degrees.
2. View TropicalCyclone intensity map.

Additional context
Add any other context or screenshots about the feature request here.

Screenshots
If applicable, add screenshots to help explain your problem.

Additional context
Add any other context about the problem here.
From the view of Tokyo, we can see the expected versus actual intensity profiles along the direction of travel of the tropical cyclone, which is the red line. Because we seem to have large time steps between cyclone footprints, we have spatial artifacts in the form of these dips instead of a smooth profile. Is there any way to get continuity here?

When we look at the larger view of Japan, we have an area with lots of artifacts within a high intensity area. Is there any way to smooth out that whole region to be roughly even intensity, since the dips in intensity are accidental and not essential.

[peanutfun]

From #916:

I can give a few of them a try and if anything useful comes out of it, what would be a good way to share back with you?

@haicyu It would be great if you could share your results here. I suggest we compare the non-smoothed to the smoothed hazard intensity, for a single track and varying time steps, visually. Would that be possible?

[haicyu]

Sure, I am currently working on another project in CLIMADA, related to River Flooding, but when I get free I will come back to this.

In the mean time I wonder if I could trouble you with a few questions from my initial investigation:

• Is there a quick way in CLIMADA to generate a 2D np array with the intensity values rasterized? Currently it seems to me that I would have to manually look up lat, lon, and intensity by index, form tuples, and then manually assign them to a 2D numpy array that represents the image we see when we use the method .plot_intensity().
• Could I use the first row of the .intensity field as the windspeed, instead of manually computing them from save windfields? Do you happen to know if they are the same? It seems to me that the plot_intensity() method is using data from the .intensity field, but I could be mistaken. The .intensity field also has five rows of length len(centroids), it seems that row [0] is corresponding to the image that we see based on a few experiments that I ran - I couldn't find any documentation on the other rows so I'm just curious if you happened to know what they represented.

[peanutfun]

From #916:

Do you have any high level considerations that would be helpful to take into account?

@haicyu For me, this relates mostly to ease of use and implementation. As I said, would be happy to incorporate a suitable smoothing algorithm. However, Climada has been stable for quite some time and the lack of such a tool was never a major issue. I would therefore make its application optional, e.g. via calling an additional method after constructing the wind field, or via an additional parameter to a current method.

I assume such a smoothing algorithm requires the full wind field at every time step and position. These are computed within the compute_windfields_sparse function, see

climada_python/climada/hazard/trop_cyclone/trop_cyclone_windfields.py

Line 992 in 46c4852

intensity = np.linalg.norm(windfields, axis=-1).max(axis=0)

The full wind field data is then discarded unless the user chooses to store_windfields. In any case, the data is then only available as a sparse matrix. What type of input data do your algorithms need? Do we need to convert the data into a dense numpy array or even an xarray DataArray?

[haicyu]

You are right @chahank that increasing time resolution can help to create smoothness in the direction of the TC's travel and thank you for your additional perspective, however as @peanutfun points out it is costly. But then from a practical perspective there remains the question of the gaps that emerge in between two parallel TC tracks, as you can see in the pictures that I have attached.

The dips in intensity, while representing a physical simulation, do not represent the climate risk that is being faced at each location. From an industry perspective, this makes it unusable. I am aware that suboptimal technical approaches with consequential tradeoffs are being taken in order to achieve a degree of smoothness.

While the workaround may be to increase the number of synthetic tracks generated, this is proving to be too costly to be effective. Focusing on just one sub-region of the world and at the limits of my computational resources, I was able to generate the attached images which still suffer from these artifacts which make it not suitable for our use case.

Smoothing is computationally far cheaper to produce the same result - by smoothing a few tracks, you could have the same effect as having simulated hundreds of tracks that still may not adequately fill those inter-track gaps, and it allows you to operate at a lower time step resolution while achieving the same result.

But I share your concern for the physicality of smoothing, particularly at the land-sea boundary, and I would be keen to learn additional concerns that you may have and brainstorm ways to address it. I will likely do this smoothing work in the coming weeks, certainly adequate from a risk management perspective and a significant improvement on the local state of the art, but I thought it might be nice to share it with others so we all have it, it feels like a nice thing to do.

[haicyu]

Oh I just saw your final comment. Sure, happy to put something together - it may be in a few weeks time but it may be more clear when we are talking in more detail. Wish you all the best in the meantime!

[chahank]

But then from a practical perspective there remains the question of the gaps that emerge in between two parallel TC tracks, as you can see in the pictures that I have attached.

The dips in intensity, while representing a physical simulation, do not represent the climate risk that is being faced at each location. From an industry perspective, this makes it unusable. I am aware that suboptimal technical approaches with consequential tradeoffs are being taken in order to achieve a degree of smoothness.

This is an interesting point to me, but I am a bit confused as to what the smoothing then could be used for. If you want to assess risk using a probabilistic event-based approach, you cannot just smooth out tracks, as this will destroy your event statistical analysis such as the average annual risk or the 1/100 year risk (besides being unphysical). What kind of results are you interested in that would benefit from a smoothing?

[chahank]

Maybe a more practical approach would be to compute a large number of tracks at lower resolution, and then use an up-scaling algorithm.

[peanutfun]

But then from a practical perspective there remains the question of the gaps that emerge in between two parallel TC tracks, as you can see in the pictures that I have attached.

@haicyu It seems I misunderstood your initial feature request. I thought you meant the gaps in the wind field of a single track, that emerge when time steps are so large that the cyclone travels several centroids within each time step. I think that these can indeed be smoothed out by an algorithm and actually improve physical consistency to some degree.

I did not realize you meant the overall coverage of the synthetic tracks. Note that generating a climatological, probabilistic TC track set is a major scientific undertaking and there are several approaches to it. The one implemented in CLIMADA is likely not competitive to the state of the art. Depending on your specific use case, it might be better to look into existing TC hazard datasets like e.g. STORM https://doi.org/10.1038/s41597-020-0381-2

[haicyu]

@peanutfun would you on the CLIMADA side have any resources available for a review? Would this be something that could be published in a paper?

[peanutfun]

What exactly do you mean by resources? We should have enough time to spare for a thorough code review. This usually happens within a GitHub pull request, see our contribution guide.

In the past, we have published submodules of Climada in journals like Geoscientific Model Development (see https://doi.org/10.5194/gmd-15-7177-2022) and we have a calibration submodule currently under review in the Journal of Open Source Software (see https://joss.theoj.org/papers/50845e31c6cb6894baae492bfd853671). The latter might be a suitable place to publish the entire TC submodule with the different wind field models implemented. However, we would have to start a conversation with the original module authors about that.

[bguillod]

To try and summarize, there are three sources of lack of smoothness in your example @haicyu, which I'm very familiar with (most of those have already been mentioned by @peanutfun or @chahank in some flavor; I'm just trying to summarize as an overview here):

1. Along a single track: If the temporal resolution of the track was relatively low before calculating the wind map. I personally typically use tracks interpolated to 10 minutes or lower (time_res_h=1/6).
2. Between each set of two tracks: this can be partly solved by generating more synthetic tracks as you suggested, however there is a limit to what this can do, especially if the spatial resolution of your centroids is very high (which by the way becomes fake precision to me when you go higher resolution than about 5km or so, because then one would need to include additional information in the wind speed modelling such as roughness length and topography - note that the "5km" mentioned here come from my intuition only).
3. However, the main issue in your example, IMHO, is that you are looking at the maximum intensity of the event set. This of course is very unstable and sensitive to individual tracks. Instead I recommend you to look at return period values instead, which are less sensitive - provided you have enough data (total number of simulated years) to robustly estimate the desired return values. For instance you could generate 1000 years of data (25 ensemble members with 40 historical years or something like this) and then estimate return values for return periods of 10, 20 or 50 years rather robustly.

I agree with @chahank that smoothing the values via an algorithm can be quite misleading in this case so I'd be cautious with such procedures.

[bguillod]

To calculate the local return value for a given return period you can use the Hazard's method local_exceedance_inten, or the plot_rp_intensity method to directly plot it.
If you e.g. select a return period of 100 (years), you get an estimate of the wind speed which would occur with a probability of 1/100 every year. This is much more robust than the maximum of all the values; moreover you can select to look at occasional events (e.g. for a return period of 10 or 20 years), rare events (e.g. 50 years) or very rare ones (e.g. 100 years or higher) (n.b. subjectivity is involved in the choice of the appropriate return period, but this is still much better than to rely on the maximum, since a maximum could be something happening every million years on average, which is probably overkill unless you're tackling something like nuclear power plant safety; also, the maximum is always the maximum of the data subset you are using and is likely to display a lower value if this data consists of 20 years compared to if it consists of 20'000 stochastically simulated years).

[haicyu]

Hi @bguillod,

I tried this method as you suggested and the results look good. Can these return period intensity maps be used with impact functions to calculate AAL?

Thanks,
Charles

[bguillod]

I tried this method as you suggested and the results look good.

Great!

Can these return period intensity maps be used with impact functions to calculate AAL?

In theory yes, you could approximate the AAL using return period maps. However, in this case it seems a convoluted way to do it: the most natural but also straightforward way to calculate the average annual loss (AAL) is to calculate the impact from your hazard data (TropCyclone) and then extract the AAL (in climada, the AAL at each exposure location is the attribute eai_exp from an Impact object, while the AAL of the whole portfolio is attribute aai_agg).

[haicyu]

to calculate the impact from your hazard data

Yes, this is what I have done but it appears to be quite sensitive to location. When I set a portfolio with a single exposure and move it around where there is a storm footprint, I get a high rate of change of AAL with respect to space, hence what brought me here with the need for smoothing. Is there any way to compute a "smoothed" AAL? That would be the end goal for me.

[bguillod]

When I set a portfolio with a single exposure and move it around where there is a storm footprint, I get a high rate of change of AAL with respect to space

Is this when using a full probabilistic event set with many years, or are you just using a few events for this? I would expect that with a large enough event set (i.e. one which is built with many years of historical tracks and enough synthetic tracks per historical track and all tracks being interpolated to a sufficient temporal resolution before calculating the wind fields) the AAL should be relatively smooth.

It is difficult to judge without looking into your data in detail, however such an investigation seems out of scope of this github issue, unless it relates to an issue in the climada code.