Calculate monthly mean data from daily inputs - differences between CDAT vs xCDAT/xarray¶

notes¶

imports¶

In [4]:

%%time
# using xcd061cdmcdu env
import cartopy.crs as ccrs
import datetime as dt
import matplotlib.pyplot as plt
import numpy as np
import os
import xarray as xr
import xcdat as xc
import cdms2 as cdm
import cdutil as cdu
import cdtime as cdt
import cftime as cft

CPU times: user 19 µs, sys: 4 µs, total: 23 µs
Wall time: 26.5 µs

set data paths¶

In [5]:

obsPath = "/p/user_pub/PCMDIobs/obs4MIPs_input/RSS/RSS-MW5-1/v20230605/"

cdat reads¶

In [6]:

print(obsPath.split("/")[1:])
obsPathXml = os.path.join("/", *obsPath.split("/")[1:7], "199801.xml")
# open file handle
fH = cdm.open(obsPathXml)
# read sst variable
sst199801 = fH("analysed_sst")
print("sst199801.getTime():", sst199801.getTime())
time199801 = sst199801.getTime().asComponentTime()
#print("time199801:", time199801)
# assign correct bounds for daily data
cdu.setTimeBoundsDaily(sst199801)
time199801d = sst199801.getTime().asComponentTime()
#print("time199801d:", time199801d)  # identical to time199801
# calculate monthly mean
sst199801mean = cdu.JAN(sst199801)
# query array shapes
print("sst199801.shape:", sst199801.shape)
print("sst199801mean.shape:", sst199801mean.shape)
# query cdat-generated time values
time199801mean = sst199801mean.getTime().asComponentTime()
print()
print("time199801mean:\n", time199801mean[0])
sst1998meanTimeBounds = sst199801mean.getTime().getBounds()[0]
# map back to relative
origin = dt.datetime(1981, 1, 1, 0, 0, 0)
startBounds = origin + dt.timedelta(0, sst1998meanTimeBounds[0])
endBounds = origin + dt.timedelta(0, sst1998meanTimeBounds[1])
print("time199801mean bounds:\n", startBounds, "\n", endBounds)
fH.close()

['p', 'user_pub', 'PCMDIobs', 'obs4MIPs_input', 'RSS', 'RSS-MW5-1', 'v20230605', '']
sst199801.getTime():    id: time
   Designated a time axis.
   units:  seconds since 1981-01-01 00:00:00
   Length: 31
   First:  536500800.0
   Last:   539092800.0
   Other axis attributes:
      axis: T
      calendar: gregorian
      realtopology: linear
   Python id:  0x7f73fb76ebc0

sst199801.shape: (31, 720, 1440)
sst199801mean.shape: (1, 720, 1440)

time199801mean:
 1998-1-16 12:0:0.0
time199801mean bounds:
 1998-01-01 00:00:00 
 1998-02-01 00:00:00

xcdat reads¶

In [7]:

def setCalendar(ds):
    # https://github.com/pydata/xarray/issues/6259
    ds.time.attrs["calendar"] = "standard"
    ds.time.attrs["units"] = "seconds since 1981-01-01 00:00:00"
    return ds
    #return xr.decode_cf(ds)

dataPath = os.path.join(obsPath, "199801*.nc")
#dataPath = os.path.join(obsPath, "1998*.nc")
print("dataPath:", dataPath)
ds = xc.open_mfdataset(dataPath, preprocess=setCalendar)
print("done!")

dataPath: /p/user_pub/PCMDIobs/obs4MIPs_input/RSS/RSS-MW5-1/v20230605/199801*.nc
done!

view dataset¶

In [8]:

ds

Out[8]:

<xarray.Dataset>
Dimensions:           (lat: 720, lon: 1440, time: 31, bnds: 2)
Coordinates:
  * lat               (lat) float32 -89.88 -89.62 -89.38 ... 89.38 89.62 89.88
  * lon               (lon) float32 -179.9 -179.6 -179.4 ... 179.4 179.6 179.9
  * time              (time) object 1998-01-01 12:00:00 ... 1998-01-31 12:00:00
Dimensions without coordinates: bnds
Data variables:
    analysed_sst      (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    analysis_error    (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    sea_ice_fraction  (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    mask              (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    lon_bnds          (lon, bnds) float32 -180.0 -179.8 -179.8 ... 179.8 180.0
    lat_bnds          (lat, bnds) float32 -90.0 -89.75 -89.75 ... 89.75 90.0
Attributes: (12/46)
    Conventions:                CF-1.8,ACDD-1.3
    title:                      Analysed foundation sea surface temperature o...
    summary:                    A merged, multi-sensor L4 foundation SST prod...
    references:                 http://www.remss.com/measurements/sea-surface...
    institution:                REMSS
    history:                    2021-11-11 18:18:15+0000 created by sst_fusio...
    ...                         ...
    project:                    Group for High Resolution Sea Surface Tempera...
    publisher_name:             The GHRSST Project Office
    publisher_email:            ghrsst-po@nceo.ac.uk
    publisher_url:              http://www.ghrsst.org
    processing_level:           L4
    cdm_data_type:              grid

add calendar attribute to time axis¶

In [9]:

ds.time.attrs['calendar'] = 'standard'

In [10]:

ds

Out[10]:

<xarray.Dataset>
Dimensions:           (lat: 720, lon: 1440, time: 31, bnds: 2)
Coordinates:
  * lat               (lat) float32 -89.88 -89.62 -89.38 ... 89.38 89.62 89.88
  * lon               (lon) float32 -179.9 -179.6 -179.4 ... 179.4 179.6 179.9
  * time              (time) object 1998-01-01 12:00:00 ... 1998-01-31 12:00:00
Dimensions without coordinates: bnds
Data variables:
    analysed_sst      (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    analysis_error    (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    sea_ice_fraction  (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    mask              (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    lon_bnds          (lon, bnds) float32 -180.0 -179.8 -179.8 ... 179.8 180.0
    lat_bnds          (lat, bnds) float32 -90.0 -89.75 -89.75 ... 89.75 90.0
Attributes: (12/46)
    Conventions:                CF-1.8,ACDD-1.3
    title:                      Analysed foundation sea surface temperature o...
    summary:                    A merged, multi-sensor L4 foundation SST prod...
    references:                 http://www.remss.com/measurements/sea-surface...
    institution:                REMSS
    history:                    2021-11-11 18:18:15+0000 created by sst_fusio...
    ...                         ...
    project:                    Group for High Resolution Sea Surface Tempera...
    publisher_name:             The GHRSST Project Office
    publisher_email:            ghrsst-po@nceo.ac.uk
    publisher_url:              http://www.ghrsst.org
    processing_level:           L4
    cdm_data_type:              grid

add missing time_bnds¶

In [11]:

ds = ds.bounds.add_missing_bounds(axes="T")
#ds.time.attrs["bounds"] = "time_bnds"

In [12]:

ds

Out[12]:

<xarray.Dataset>
Dimensions:           (lat: 720, lon: 1440, time: 31, bnds: 2)
Coordinates:
  * lat               (lat) float32 -89.88 -89.62 -89.38 ... 89.38 89.62 89.88
  * lon               (lon) float32 -179.9 -179.6 -179.4 ... 179.4 179.6 179.9
  * time              (time) object 1998-01-01 12:00:00 ... 1998-01-31 12:00:00
Dimensions without coordinates: bnds
Data variables:
    analysed_sst      (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    analysis_error    (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    sea_ice_fraction  (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    mask              (time, lat, lon) float32 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    lon_bnds          (lon, bnds) float32 -180.0 -179.8 -179.8 ... 179.8 180.0
    lat_bnds          (lat, bnds) float32 -90.0 -89.75 -89.75 ... 89.75 90.0
    time_bnds         (time, bnds) object 1998-01-01 00:00:00 ... 1998-02-01 ...
Attributes: (12/46)
    Conventions:                CF-1.8,ACDD-1.3
    title:                      Analysed foundation sea surface temperature o...
    summary:                    A merged, multi-sensor L4 foundation SST prod...
    references:                 http://www.remss.com/measurements/sea-surface...
    institution:                REMSS
    history:                    2021-11-11 18:18:15+0000 created by sst_fusio...
    ...                         ...
    project:                    Group for High Resolution Sea Surface Tempera...
    publisher_name:             The GHRSST Project Office
    publisher_email:            ghrsst-po@nceo.ac.uk
    publisher_url:              http://www.ghrsst.org
    processing_level:           L4
    cdm_data_type:              grid

try xCDAT temporal.group_average¶

In [13]:

mean_monthXc = ds.temporal.group_average("analysed_sst", freq="month", weighted=True)

In [14]:

mean_monthXc
# loses bounds
# time should be 1998-1-16 12:00, but is rather 1998-1-1 00:00

Out[14]:

<xarray.Dataset>
Dimensions:       (lat: 720, lon: 1440, bnds: 2, time: 1)
Coordinates:
  * lat           (lat) float32 -89.88 -89.62 -89.38 ... 89.38 89.62 89.88
  * lon           (lon) float32 -179.9 -179.6 -179.4 ... 179.4 179.6 179.9
  * time          (time) object 1998-01-01 00:00:00
Dimensions without coordinates: bnds
Data variables:
    lon_bnds      (lon, bnds) float32 -180.0 -179.8 -179.8 ... 179.8 179.8 180.0
    lat_bnds      (lat, bnds) float32 -90.0 -89.75 -89.75 ... 89.75 89.75 90.0
    analysed_sst  (time, lat, lon) float64 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
Attributes: (12/46)
    Conventions:                CF-1.8,ACDD-1.3
    title:                      Analysed foundation sea surface temperature o...
    summary:                    A merged, multi-sensor L4 foundation SST prod...
    references:                 http://www.remss.com/measurements/sea-surface...
    institution:                REMSS
    history:                    2021-11-11 18:18:15+0000 created by sst_fusio...
    ...                         ...
    project:                    Group for High Resolution Sea Surface Tempera...
    publisher_name:             The GHRSST Project Office
    publisher_email:            ghrsst-po@nceo.ac.uk
    publisher_url:              http://www.ghrsst.org
    processing_level:           L4
    cdm_data_type:              grid

In [15]:

#mean_monthXc2 = mean_monthXc.bounds.add_time_bounds(method="freq", freq="month")
# reassigning bounds fails due to 
# ValueError: Cannot generate bounds for coordinate variable 'time' which has a length <= 1 (singleton).

correct xCDAT temporal.group_average¶

In [16]:

#mean_monthXc.time
# update new time value
mean_monthXc.variables["time"]
newTime = [cft.DatetimeGregorian(1998, 1, 16, 12, 0, 0, 0, has_year_zero=False)]
mean_monthXc.update({"time": ("time", newTime)})
# assign new time_bnds variable
newTimeBoundVals = np.array([[cft.DatetimeGregorian(1998, 1, 1, 0, 0, 0, 0, has_year_zero=False),
                 cft.DatetimeGregorian(1998, 2, 1, 0, 0, 0, 0, has_year_zero=False)]])
print(type(newTimeBoundVals))
print(newTimeBoundVals.shape)
newTimeBounds = xr.Dataset({"time_bnds": (("time", "bnds"), newTimeBoundVals)})
#newTimeBounds
mean_monthXc.assign(newTimeBounds)

<class 'numpy.ndarray'>
(1, 2)

Out[16]:

<xarray.Dataset>
Dimensions:       (lat: 720, lon: 1440, bnds: 2, time: 1)
Coordinates:
  * lat           (lat) float32 -89.88 -89.62 -89.38 ... 89.38 89.62 89.88
  * lon           (lon) float32 -179.9 -179.6 -179.4 ... 179.4 179.6 179.9
  * time          (time) object 1998-01-16 12:00:00
Dimensions without coordinates: bnds
Data variables:
    lon_bnds      (lon, bnds) float32 -180.0 -179.8 -179.8 ... 179.8 179.8 180.0
    lat_bnds      (lat, bnds) float32 -90.0 -89.75 -89.75 ... 89.75 89.75 90.0
    analysed_sst  (time, lat, lon) float64 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    time_bnds     (time, bnds) object 1998-01-01 00:00:00 1998-02-01 00:00:00
Attributes: (12/46)
    Conventions:                CF-1.8,ACDD-1.3
    title:                      Analysed foundation sea surface temperature o...
    summary:                    A merged, multi-sensor L4 foundation SST prod...
    references:                 http://www.remss.com/measurements/sea-surface...
    institution:                REMSS
    history:                    2021-11-11 18:18:15+0000 created by sst_fusio...
    ...                         ...
    project:                    Group for High Resolution Sea Surface Tempera...
    publisher_name:             The GHRSST Project Office
    publisher_email:            ghrsst-po@nceo.ac.uk
    publisher_url:              http://www.ghrsst.org
    processing_level:           L4
    cdm_data_type:              grid

try xCDAT temporal.average¶

In [17]:

mean_monthXc2 = ds.temporal.average("analysed_sst", weighted=True)

In [18]:

mean_monthXc2
# loses time dimension altogether

Out[18]:

<xarray.Dataset>
Dimensions:       (lat: 720, lon: 1440, bnds: 2)
Coordinates:
  * lat           (lat) float32 -89.88 -89.62 -89.38 ... 89.38 89.62 89.88
  * lon           (lon) float32 -179.9 -179.6 -179.4 ... 179.4 179.6 179.9
Dimensions without coordinates: bnds
Data variables:
    lon_bnds      (lon, bnds) float32 -180.0 -179.8 -179.8 ... 179.8 179.8 180.0
    lat_bnds      (lat, bnds) float32 -90.0 -89.75 -89.75 ... 89.75 89.75 90.0
    analysed_sst  (lat, lon) float64 dask.array<chunksize=(720, 1440), meta=np.ndarray>
Attributes: (12/46)
    Conventions:                CF-1.8,ACDD-1.3
    title:                      Analysed foundation sea surface temperature o...
    summary:                    A merged, multi-sensor L4 foundation SST prod...
    references:                 http://www.remss.com/measurements/sea-surface...
    institution:                REMSS
    history:                    2021-11-11 18:18:15+0000 created by sst_fusio...
    ...                         ...
    project:                    Group for High Resolution Sea Surface Tempera...
    publisher_name:             The GHRSST Project Office
    publisher_email:            ghrsst-po@nceo.ac.uk
    publisher_url:              http://www.ghrsst.org
    processing_level:           L4
    cdm_data_type:              grid

xarray reads¶

In [19]:

mean_monthXr = ds.analysed_sst.resample(time='ME').mean()
# https://stackoverflow.com/questions/50564459/using-xarray-to-make-monthly-average

In [20]:

mean_monthXr
# time should be 1998-1-16 12:00, but is rather 1998-1-31 00:00

Out[20]:

<xarray.DataArray 'analysed_sst' (time: 1, lat: 720, lon: 1440)>
dask.array<stack, shape=(1, 720, 1440), dtype=float32, chunksize=(1, 720, 1440), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float32 -89.88 -89.62 -89.38 -89.12 ... 89.38 89.62 89.88
  * lon      (lon) float32 -179.9 -179.6 -179.4 -179.1 ... 179.4 179.6 179.9
  * time     (time) object 1998-01-31 00:00:00
Attributes:
    units:                  K
    long_name:              analysed sea surface temperature
    standard_name:          sea_surface_foundation_temperature
    valid_min:              -32767
    valid_max:              32767
    source:                 REMSS-L3C-TMI
    coverage_content_type:  physicalMeasurement

In [21]:

ds.analysed_sst

Out[21]:

<xarray.DataArray 'analysed_sst' (time: 31, lat: 720, lon: 1440)>
dask.array<concatenate, shape=(31, 720, 1440), dtype=float32, chunksize=(1, 720, 1440), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float32 -89.88 -89.62 -89.38 -89.12 ... 89.38 89.62 89.88
  * lon      (lon) float32 -179.9 -179.6 -179.4 -179.1 ... 179.4 179.6 179.9
  * time     (time) object 1998-01-01 12:00:00 ... 1998-01-31 12:00:00
Attributes:
    units:                  K
    long_name:              analysed sea surface temperature
    standard_name:          sea_surface_foundation_temperature
    valid_min:              -32767
    valid_max:              32767
    source:                 REMSS-L3C-TMI
    coverage_content_type:  physicalMeasurement

In [22]:

xr.show_versions()

INSTALLED VERSIONS
------------------
commit: None
python: 3.10.13 | packaged by conda-forge | (main, Dec 23 2023, 15:36:39) [GCC 12.3.0]
python-bits: 64
OS: Linux
OS-release: 3.10.0-1160.90.1.el7.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: en_US.UTF-8
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.14.3
libnetcdf: 4.9.2

xarray: 2023.12.0
pandas: 2.1.4
numpy: 1.26.3
scipy: 1.11.4
netCDF4: 1.6.5
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.6.3
nc_time_axis: None
iris: None
bottleneck: 1.3.7
dask: 2023.12.1
distributed: 2023.12.1
matplotlib: 3.8.2
cartopy: 0.22.0
seaborn: None
numbagg: None
fsspec: 2023.12.2
cupy: None
pint: None
sparse: 0.15.1
flox: None
numpy_groupies: None
setuptools: 69.0.3
pip: 23.3.2
conda: None
pytest: None
mypy: None
IPython: 8.20.0
sphinx: None

/home/durack1/mambaforge/envs/xcd061cdmcdu/lib/python3.10/site-packages/_distutils_hack/__init__.py:33: UserWarning: Setuptools is replacing distutils.
  warnings.warn("Setuptools is replacing distutils.")

In [ ]: