A python library to evaluate the skill of n-dimension data
There are lots of python libraries out there which provide metrics of skill. Metrics of skill are error calculations which tells you how good or how bad a model was compared to what actually happened. scikit-learn and scipy offer a variety of metrics as can be found here https://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics and https://docs.scipy.org/doc/scipy/reference/stats.html. However, most of them work with one-dimensional data: data which changes with time, for example.
In geosciences we often deal with n-dimensional data. For example, atmospheric data is traditionally four dimensional: latitude, longitude, height and time. An atmospheric model is often assessed using the time dimension. How does the forecast error grow over time? Therefore, the model data for the dimensions: latitude, longitude and height are being compared to the observed data with latitude, longitude and height, over time. However, for other data you may want to assess skill over another dimension. I wanted to capture this flexibility in the code. I also wanted to assure the calculations remain vectorized to speed up computation. Looping over latitude then looping over longitude then looping over height and applying the one dimension skill metrics above would be computationally expensive.
Xarray is a fantastic python library which labels n-dimensional data and allow for you to easily perform operations on a particular dimension. If you data has the four dimensions stated above and you wanted to do what was the maximum over a particular dimension you would do:
da.max(dim='time)Where da is an Xarray.DataArray and contains the data and labels. This will reduce the label down to three dimensions: latitude, longitude and height.
Xarray has an apply_ufunc method to extend the capability of xarray. This allows for more specific xarray functions that are too broad to be included in the package. It is here where you write a function for numpy arrays (similar numeric array) and then extend it to Xarray data objects. So xskillscore has two parts. First, I write the skill metric for n-dimensional numpy arrays. The key for this was to use numpy.rollaxis to put the axis I am assessing skill over (e.g. time) first. You can take a look at the code here: https://github.com/raybellwaves/xskillscore/blob/master/xskillscore/core/np_deterministic.py. Secondly, you put this function into xarray.apply_ufunc so it will be applied to two xarray.DataArray's (the model data and the observed data).
You can use the package as
r = xs.pearson_r(da_obs, da_fct, 'time')if da_obs and da_fct and four dimensional as stated above, r will be the correlation coefficient with shape latitude, longitude and height. Similarly you could do
r = xs.pearson_r(da_obs, da_fct, 'height')which will return a three dimension xarray.DataArray with dimensions latitude, longitude and time.
I only got up to including Pearson's correlation coefficient, the 2-tailed p-value associated with pearson's correlation coefficient and the Root Mean Squared Error in the package. These are deterministic metrics. That is a metric which do not take into account the distribution of the forecast such as if you have ensemble information. If your forecasts contained ensembles you would use the mean of those forecasts. However, it was all that was required to produce a figure in a BAMS paper which recently got submitted on the SubX project.
A big thank you to the xarray developers for the package and there guidance in helping me put this together (https://groups.google.com/forum/#!topic/xarray/z8ue0G-BLc8). I'm happy this package is listed on the xarray docs http://xarray.pydata.org/en/stable/related-projects.html#extend-xarray-capabilities.