CDAT / cdms

8 stars 10 forks source link

MV2.repeat unexpected behavior #423

Open pochedls opened 3 years ago

pochedls commented 3 years ago

Describe the bug Behavior of MV2.repeat is unexpected. If I have a transientVariable a = [0, 1, 2, 3], and apply MV2.repeat(a, 2) I get: [0, 0, 1, 1, 2, 2, 3, 3] whereas I would expect [0, 1, 2, 3, 0, 1, 2, 3].

Is MV2.repeat doing the right thing?

To Reproduce

Simple example

[In]: import cdms2
[In]: import MV2
[In]: a = cdms2.createVariable([0, 1, 2, 3])
[In]: MV2.repeat(a, 2)
[Out]: masked_array(data=[0, 0, 1, 1, 2, 2, 3, 3],
             mask=False,
       fill_value=999999)

Applied Example: Remove seasonal cycle from each gridpoint

import cdms2
import MV2
import numpy as np
import matplotlib.pyplot as plt

fn = '/p/user_pub/work/pochedley1/ssmi/ssmi_vapor_1Deg.nc3.nc'
f = cdms2.open(fn)
prw = f('prw')
f.close()
time = prw.getTime()
decimalTime = np.arange(1988, 2020, 1/12.)
lat = prw.getLatitude()
lon = prw.getLongitude()
nyrs = int(len(time)/12)
tropicalInds = np.where(np.abs(lat) < 20)[0]

# MV2 to subtract seasonal cycle
prwClimatology = MV2.reshape(prw, (nyrs, 12, len(lat), len(lon)))  # reshape to (year, month, lat, lon)
prwClimatology = MV2.average(prwClimatology, axis=0)  # average over all years
prwAnom = prw - MV2.repeat(prwClimatology, nyrs, axis=0)  # substract repeated climatology
# simple tropical average (for plotting)
prwAnomTrop = MV2.average(MV2.average(MV2.take(prwAnom, tropicalInds, axis=1), axis=2), axis=1)

# numpy to subtract seasonal cycle
prw_n = np.array(prw)  # to numpy array
prw_n[prw_n < 0] = np.nan  # nan missing data
prw_n[prw_n > 1000] = np.nan  # nan missing data
prwClimatology_n = np.reshape(prw_n, (nyrs, 12, len(lat), len(lon)))  # reshape to (year, month, lat, lon)
prwClimatology_n = np.mean(prwClimatology_n, axis=0)  # average over all years
prwAnom_n = prw_n - np.tile(prwClimatology_n, (nyrs, 1, 1))  # substract repeated climatology
# simple tropical average (for plotting)
prwAnomTrop_n = np.nanmean(np.nanmean(np.take(prwAnom_n, tropicalInds, axis=1), axis=2), axis=1)

# plot two series
plt.plot(decimalTime, np.array(prwAnomTrop), label='MV2')
plt.plot(decimalTime, prwAnomTrop_n, label='numpy')
plt.legend()
plt.xlabel('Time')
plt.ylabel('tropical prw')
plt.show()

test

Note: the numpy version of the time series is correct and the expected behavior.

Expected behavior In the longer example, I expect MV2.repeat to behave like np.tile (repeating the annual cycle so that it can be subtracted from the prw time series in order to produce the anomaly time series).

Desktop (please complete the following information):