Trends with the SAGE II / OSIRIS / OMPS-LP Dataset

Here we calculate trends using the SAGE II / OSIRIS / OMPS-LP dataset by regressing to deseasonalized relative monthly zonal mean anomalies

import xarray as xr
import numpy as np
from LOTUS_regression.regression import regress_all_bins
from LOTUS_regression.predictors import load_data

The data is in a single NetCDF4 file.

MERGED_FILE = r'/home/runner/work/lotus-regression/lotus-regression/test_data//S2_OS_OMPS/MERGED_LOTUS.nc'

mzm_data = xr.open_dataset(MERGED_FILE, engine='netcdf4')

print(mzm_data)

<xarray.Dataset> Size: 14MB
Dimensions:           (altitude: 59, time: 387, mean_latitude: 13, month: 12)
Coordinates:
  * altitude          (altitude) float64 472B 0.0 1.0 2.0 3.0 ... 56.0 57.0 58.0
  * time              (time) datetime64[ns] 3kB 1984-10-01 ... 2016-12-01
    latitude_bin      (mean_latitude) <U14 728B ...
  * mean_latitude     (mean_latitude) float64 104B -60.0 -50.0 ... 50.0 60.0
  * month             (month) int64 96B 1 2 3 4 5 6 7 8 9 10 11 12
Data variables:
    anomaly           (mean_latitude, time, altitude) float64 2MB ...
    relative_anomaly  (mean_latitude, time, altitude) float64 2MB ...
    std               (mean_latitude, time, altitude) float64 2MB ...
    relative_std      (mean_latitude, time, altitude) float64 2MB ...
    count             (mean_latitude, time, altitude) float64 2MB ...
    osiris_mzm        (mean_latitude, month, altitude) float64 74kB ...
    sageii_mzm        (mean_latitude, month, altitude) float64 74kB ...
    omps_mzm          (mean_latitude, month, altitude) float64 74kB ...
    inst_flags        (mean_latitude, time, altitude) float64 2MB ...

predictors = load_data('pred_baseline_pwlt.csv')

print(predictors.columns)

Index(['enso', 'solar', 'qboA', 'qboB', 'aod', 'linear_pre', 'linear_post',
       'constant'],
      dtype='object')

mzm_data = mzm_data.sel(mean_latitude=slice(-60, 60), altitude=slice(10, 50))  # remove bins without data
results = regress_all_bins(predictors, mzm_data['relative_anomaly'], tolerance=0.1)

# Convert to ~ percent
results *= 100

print(results)

<xarray.Dataset> Size: 69kB
Dimensions:          (mean_latitude: 13, altitude: 41)
Coordinates:
  * mean_latitude    (mean_latitude) float64 104B -60.0 -50.0 ... 50.0 60.0
  * altitude         (altitude) float64 328B 10.0 11.0 12.0 ... 48.0 49.0 50.0
Data variables: (12/16)
    enso             (mean_latitude, altitude) float64 4kB -1.528 ... 0.3204
    enso_std         (mean_latitude, altitude) float64 4kB 1.173 ... 0.3949
    solar            (mean_latitude, altitude) float64 4kB -0.9302 ... 0.4314
    solar_std        (mean_latitude, altitude) float64 4kB 1.215 ... 0.4279
    qboA             (mean_latitude, altitude) float64 4kB -3.226 ... -0.5349
    qboA_std         (mean_latitude, altitude) float64 4kB 1.137 ... 0.3806
    ...               ...
    linear_pre       (mean_latitude, altitude) float64 4kB 7.944 ... -6.499
    linear_pre_std   (mean_latitude, altitude) float64 4kB 4.169 3.295 ... 1.361
    linear_post      (mean_latitude, altitude) float64 4kB -12.62 ... 3.515
    linear_post_std  (mean_latitude, altitude) float64 4kB 2.199 1.704 ... 0.717
    constant         (mean_latitude, altitude) float64 4kB 14.9 9.552 ... -4.378
    constant_std     (mean_latitude, altitude) float64 4kB 2.684 ... 0.7894

import LOTUS_regression.plotting.trends as trends
trends.pre_post_with_confidence(results, x='mean_latitude', y='altitude', ylim=(18, 50), log_y=False, figsize=(16, 6),
                                x_label='Latitude [$^\circ$]', y_label='Altitude [km]', pre_title='Pre 1997',
                                post_title='Post 1997')