Hi
We can use the function climate.anomaly(...) to calculate anomalies:
data_hist_clim = ct.climate.climatology_mean(data_hist) data_anomalies = ct.climate.anomaly(data_scen, climatology=data_hist_clim)
From looking at the source code I gather it calculates the anomalies by subtracting a reference from the data.
But sometimes (e.g., in case of precipitation) relative anomalies (i.e. a ratio), or anomalies as a percentage of the original, make more sense. What is the best way of doing this? I can of course calculate these myself:
# first mask out very small values in the reference data data_hist_clim = ct.cube.where(data_hist_clim > 0.001, data_hist_clim) data_anomalies = 100*(data_scen - data_hist_clim) / data_hist_clim
Although this seems to work, and appears to give me sensible values, and I can even visualise the result offline, my workflow then fails when I try to visualise these data, e.g. in a LiveMap. Presumably this is because climate.anomaly(...) does something with the dimensions or attributes of the result (data_anomalies) that my simple calculation doesn't, but exactly what I cannot figure out.
So, what is the best way to calculate relative anomalies?
Many thanks!
Rutger
9 Comments
James Varndell
Hi Rutger,
Thank you for your question.
I think the method you're using for calculating anomalies as a percentage is the best way to do it at the moment. Could you give some more information about the error you get when trying to visualise these data? If you have an example workflow I could play with that would be ideal!
Many thanks,
James
Rutger Dankers
Hi James
The error I is copied below. Although the actual workflow I am working on is more complicated than this, this workflow reproduces the issue: https://cds.climate.copernicus.eu/toolbox-editor/13784/precip_anomalies_example The calculation of anomalies is done in the function calc_anom(...). Literally if I switch the calculation from using climate.anomaly(...) to percentages as outlined above, the script fails in the visualisation step.
Error message:
Traceback (most recent call last): File "/opt/cdstoolbox/cdscompute/cdscompute/cdshandlers/services/handler.py", line 49, in handle_request result = cached(context.method, proc, context, context.args, context.kwargs) File "/opt/cdstoolbox/cdscompute/cdscompute/caching.py", line 108, in cached result = proc(context, *context.args, **context.kwargs) File "/opt/cdstoolbox/cdscompute/cdscompute/services.py", line 118, in __call__ return p(*args, **kwargs) File "/opt/cdstoolbox/cdscompute/cdscompute/services.py", line 59, in __call__ return self.proc(context, *args, **kwargs) File "/home/cds/cdsservices/services/wms/livemap_legacy.py", line 440, in plot land_style, coastline_style) File "/home/cds/cdsservices/services/wms/livemap.py", line 196, in plot crs=crs, **kwargs) File "/opt/cdstoolbox/cdscompute/cdscompute/download.py", line 154, in to_wms self._wms = self._to_wms(*args, **kwargs) File "/opt/cdstoolbox/cdscompute/cdscompute/download.py", line 221, in _to_wms engine=self._engine, **kwargs) File "/opt/cdstoolbox/cdscdm/cdscdm/wms/layers.py", line 335, in from_wms_related_data return _animated_geojson_from_data(context, dataarray, shape, style=style, **kwargs) File "/opt/cdstoolbox/cdscdm/cdscdm/wms/layers.py", line 410, in _animated_geojson_from_data dataarray.sel(**{dim_label: fid}).values.astype(float)] File "/opt/cdstoolbox/cdscdm/cdscdm/wms/layers.py", line 409, in <listcomp> values = [i if np.isfinite(i) else None for i in ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()Thanks!
Rutger
James Varndell
Hi Rutger,
It looks like you've got both a
timeand amonthdimension in your data, which is confusing the livemap (sorry for the poor error message - this needs to be caught more usefully!).I've 'fixed' your workflow by selecting just the first
monthindex, but I'll leave it up to you how best to handle themonthdimension. Here's a working version of your workflow (the only line I've added are lines 145 and 146): https://cds.climate.copernicus.eu/toolbox-editor/20413/precip_anomalies_example.Hope that helps!
James
Rutger Dankers
Hi James
Thanks. But this is exactly my point, why is the extra
monthdimension problematic when I do the calculation myself, but not when I use the standard functionclimate.anomaly(...)?As far as I can see, the
monthdimension is created when calculating the climatology over the reference period. This is how my arrays look like (using the names in my example script):##### data_scen <xarray.DataArray 'tp' (time: 348, lat: 74, lon: 140)> dask.array<xarray-tp, shape=(348, 74, 140), dtype=float32, chunksize=(48, 74, 140), chunktype=numpy.ndarray> Coordinates: * time (time) datetime64[ns] 2071-01-01 2071-02-01 ... 2099-12-01 * lat (lat) float64 34.25 34.75 35.25 35.75 ... 69.25 69.75 70.25 70.75 * lon (lon) float64 -24.75 -24.25 -23.75 -23.25 ... 43.75 44.25 44.75 ##### data_hist <xarray.DataArray 'tp' (time: 360, lat: 74, lon: 140)> dask.array<xarray-tp, shape=(360, 74, 140), dtype=float32, chunksize=(48, 74, 140), chunktype=numpy.ndarray> Coordinates: * time (time) datetime64[ns] 1981-01-01 1981-02-01 ... 2010-12-01 * lat (lat) float64 34.25 34.75 35.25 35.75 ... 69.25 69.75 70.25 70.75 * lon (lon) float64 -24.75 -24.25 -23.75 -23.25 ... 43.75 44.25 44.75 ##### data_hist_clim <xarray.DataArray 'tp_climatology_mean' (month: 12, lat: 74, lon: 140)> [124320 values with dtype=float32] Coordinates: * lat (lat) float64 34.25 34.75 35.25 35.75 ... 69.25 69.75 70.25 70.75 * lon (lon) float64 -24.75 -24.25 -23.75 -23.25 ... 43.75 44.25 44.75 * month (month) int64 1 2 3 4 5 6 7 8 9 10 11 12 ##### data_anomalies # calculated with climate.anomaly(...) <xarray.DataArray 'tp_anomaly' (time: 348, lat: 74, lon: 140)> dask.array<xarray-tp_anomaly, shape=(348, 74, 140), dtype=float32, chunksize=(48, 74, 140), chunktype=numpy.ndarray> Coordinates: * lat (lat) float64 34.25 34.75 35.25 35.75 ... 69.25 69.75 70.25 70.75 * lon (lon) float64 -24.75 -24.25 -23.75 -23.25 ... 43.75 44.25 44.75 * time (time) datetime64[ns] 2071-01-01 2071-02-01 ... 2099-12-01 month (time) int64 dask.array<chunksize=(48,), meta=np.ndarray> Attributes: long_name: Precipitation sum anomaly units: m cell_methods: time: sum standard_name: lwe_thickness_of_precipitation_amount_anomaly type: real ##### data_anomalies # calculated manually as percentage <xarray.DataArray 'data' (time: 348, lat: 74, lon: 140, month: 12)> dask.array<xarray-data, shape=(348, 74, 140, 12), dtype=float32, chunksize=(48, 74, 140, 12), chunktype=numpy.ndarray> Coordinates: * time (time) datetime64[ns] 2071-01-01 2071-02-01 ... 2099-12-01 * lat (lat) float64 34.25 34.75 35.25 35.75 ... 69.25 69.75 70.25 70.75 * lon (lon) float64 -24.75 -24.25 -23.75 -23.25 ... 43.75 44.25 44.75 * month (month) int64 1 2 3 4 5 6 7 8 9 10 11 12 Attributes: units: 1So what exactly does
climate.anomaly(...)do to get a non-starredmonthdimension, that my manual calculation of percentages doesn't do?Thanks
Rutger
James Varndell
Hi Rutger Dankers,
Here's a quick breakdown of what's happening:
ct.climate.climatology_mean()takes some data with atimedimension and produces a climatology with amonthdimension;ct.climate.anomaly()takes some data with atimedimension and a climatology with just amonthdimension, and uses some handy xarray functionality behind the scenes to compare the monthly climatology to each month in your data in order to calculate an anomaly;timedimension and your climatology just has amonthdimension - the dimensionality of your data needs to match to do proper arithmetic.So while I think it's possible to put together a messy work around for this issue within your workflow, I think what we really need is a
relativeorpercentageargument toct.climate.anomaly, that simply produces anomalies as percentages straight out of the service. I've raised a ticket to add this argument and it should be included in next week's release (11th Nov) - does that timescale work for you? If you need something quicker, I'm happy to put together a temporary work around.Apologies that my earlier advice was wrong on this!
Many thanks,
James
Rutger Dankers
Hi James Varndell
Many thanks for your response. I'm happy to wait for a proper solution, if it can indeed be implemented by next week.
Cheers
Rutger
James Varndell
Hi Rutger,
I've just realised that there is a way to do this in the Toolbox at the moment that's not too fiddly of a workaround - you can use the
ct.cube.groupby_operatorandct.cube.groupby_reduceservices to directly compare single month values to multi-year monthly time dimensions. I've had a go at implementing this in this workflow: https://cds.climate.copernicus.eu/toolbox-editor/20413/precip_anomalies_example-1 - it seems to work, but the anomaly values are very high - is that to be expected? I might've done something slightly wrong in my calculation but hopefully this example demonstrates how you might be able to make use of thesegroupbyservices.In the meantime, we will still implement a dedicated service for calculating percentage anomalies in the Toolbox.
Thanks again,
James
Rutger Dankers
Great to see there's now a function relative_anomaly(...) to calculate relative anomalies:
James Varndell what are the units coming out of this?
Thanks
Rutger
James Varndell
Hi Rutger Dankers,
The return value from
ct.climate.relative_anomalyis the (unitless) ratio of the anomaly to the absolute value - so for example, an increase of 50% will be returned as 0.5. You can convert these values to percentages by simply multiplying the result by 100 in the Toolbox, i.e:Hope that helps!
Kind regards,
James