Skip to content

Commit

Permalink
Update climatologies to use start and end times
Browse files Browse the repository at this point in the history 
There is no longer an advantage to computing monthly climatologies
first, then seasonal climatologies (as was done in several analysis
tasks), so seasonal climatologies are now computed directly.

The indexNino34 task has been updated to work with an xarray.Dataset
instead of an xarray.DataArray because climatology calculations now
only work with the former.

The MOC time average is now computed using compute_climatology
(rather than just .mean('Time')) so that the number of days in a
month is properly accounted for.
  • Loading branch information
@xylar
xylar committed Mar 29, 2017
1 parent ab92fff commit fefb7aa
Show file tree
Hide file tree
Showing 5 changed files with 220 additions and 118 deletions.
3 changes: 2 additions & 1 deletion mpas_analysis/ocean/meridional_overturning_circulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -303,7 +303,8 @@ def _compute_moc_climo_postprocess(config, runStreams, variableMap, calendar,
dictClimo['endYearClimo'] = endYear

# Compute annual climatology
annualClimatology = ds.mean('Time')
annualClimatology = \
climatology.compute_climatology(ds, np.arange(1,13), calendar)

# Convert to numpy arrays
horizontalVel = annualClimatology.avgNormalVelocity.values
Expand Down
30 changes: 17 additions & 13 deletions mpas_analysis/ocean/nino34_index.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,9 @@
Last Modified
-------------
03/23/2017
03/29/2017
"""

import xarray as xr
import pandas as pd
import numpy as np
Expand Down Expand Up @@ -51,7 +52,7 @@ def nino34_index(config, streamMap=None, variableMap=None): # {{{
Last Modified
-------------
03/23/2017
03/29/2017
"""

print ' Load SST data...'
Expand Down Expand Up @@ -87,9 +88,8 @@ def nino34_index(config, streamMap=None, variableMap=None): # {{{
startDate=startDate,
endDate=endDate)

SSTregions = ds.avgSurfaceTemperature
print ' Compute NINO3.4 index...'
nino34 = compute_nino34_index(SSTregions[:, regionIndex], calendar)
nino34 = compute_nino34_index(ds.isel(nOceanRegions=regionIndex), calendar)

print ' Computing NINO3.4 power spectra...'
f, spectra, conf99, conf95, redNoise = compute_nino34_spectra(nino34)
Expand All @@ -108,7 +108,7 @@ def nino34_index(config, streamMap=None, variableMap=None): # {{{
# }}}


def compute_nino34_index(regionSST, calendar): # {{{
def compute_nino34_index(ds, calendar): # {{{
"""
Computes nino34 index time series. It follow the standard nino34
algorithm, i.e.,
Expand All @@ -122,15 +122,15 @@ def compute_nino34_index(regionSST, calendar): # {{{
Parameters
----------
regionSST : xarray.dataArray object
ds : xarray.DataArray object
values of SST in the nino region
calendar: {'gregorian', 'gregorian_noleap'}
The name of the calendars used in the MPAS run
Returns
-------
xarray.DataArray object containing the nino34index
xarray.Dataset object containing the nino34index
Author
------
Expand All @@ -141,15 +141,19 @@ def compute_nino34_index(regionSST, calendar): # {{{
03/23/2017
"""

if not isinstance(regionSST, xr.core.dataarray.DataArray):
raise ValueError('regionSST should be an xarray DataArray')
if not isinstance(ds, xr.core.dataset.Dataset):
raise ValueError('ds should be an xarray Dataset')

ds = climatology.add_months_and_days_in_month(ds, calendar)

# Compute monthly average and anomaly of climatology of SST
monthlyClimatology = climatology.compute_monthly_climatology(regionSST,
calendar)
anomalySST = regionSST.groupby('month') - monthlyClimatology
monthlyClimatology = climatology.compute_monthly_climatology(ds)

print monthlyClimatology

anomaly = ds.groupby('month') - monthlyClimatology

return _running_mean(anomalySST.to_pandas()) # }}}
return _running_mean(anomaly.avgSurfaceTemperature.to_pandas()) # }}}


def compute_nino34_spectra(nino34Index): # {{{
Expand Down
53 changes: 24 additions & 29 deletions mpas_analysis/ocean/ocean_modelvsobs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
Last Modified
-------------
03/23/2017
03/29/2017
"""

import xarray as xr
Expand All @@ -29,13 +29,14 @@
from ..shared.generalized_reader.generalized_reader \
import open_multifile_dataset

from ..shared.timekeeping.utility import get_simulation_start_time, \
days_to_datetime
from ..shared.timekeeping.utility import get_simulation_start_time

from ..shared.climatology import climatology

from ..shared.analysis_task import setup_task

from ..shared.mpas_xarray import mpas_xarray


def ocn_modelvsobs(config, field):

Expand All @@ -56,7 +57,7 @@ def ocn_modelvsobs(config, field):
Last Modified
-------------
03/23/2017
03/29/2017
"""

# perform common setup for the task
Expand Down Expand Up @@ -123,6 +124,8 @@ def ocn_modelvsobs(config, field):

iselvals = None

obsFieldName = 'mld_dt_mean'

if buildObsClimatologies:
# Load MLD observational data
dsObs = xr.open_mfdataset(obsFileName)
Expand All @@ -133,17 +136,18 @@ def ocn_modelvsobs(config, field):
# Rename the dimensions to be consistent with other obs. data sets
dsObs.rename({'month': 'calmonth', 'lat': 'latCoord',
'lon': 'lonCoord'}, inplace=True)
dsObs.rename({'iMONTH': 'month', 'iLAT': 'lat', 'iLON': 'lon'},
dsObs.rename({'iMONTH': 'Time', 'iLAT': 'lat', 'iLON': 'lon'},
inplace=True)
# set the coordinates now that the dimensions have the same names
dsObs.coords['lat'] = dsObs['latCoord']
dsObs.coords['lon'] = dsObs['lonCoord']
dsObs.coords['month'] = dsObs['calmonth']
dsObs.coords['Time'] = dsObs['calmonth']
dsObs['month'] = ('Time', np.array(dsObs['calmonth'], int))

dsObs = mpas_xarray.subset_variables(dsObs, [obsFieldName,
'month'])
# Reorder dataset for consistence with other obs. data sets
dsObs = dsObs.transpose('month', 'lat', 'lon')

obsFieldName = 'mld_dt_mean'
dsObs = dsObs.transpose('Time', 'lat', 'lon')

# Set appropriate MLD figure labels
observationTitleLabel = \
Expand All @@ -169,11 +173,10 @@ def ocn_modelvsobs(config, field):

if buildObsClimatologies:
dsObs = xr.open_mfdataset(obsFileName)

dsTimeSlice = dsObs.sel(time=slice(timeStart, timeEnd))
monthlyClimatology = dsTimeSlice.groupby('time.month').mean('time')

dsObs = monthlyClimatology.transpose('month', 'lat', 'lon')
dsObs.rename({'time': 'Time'}, inplace=True)
dsObs = dsObs.transpose('Time', 'lat', 'lon')
dsObs = dsObs.sel(Time=slice(timeStart, timeEnd))
dsObs['month'] = dsObs['Time.month']

obsFieldName = 'SST'

Expand All @@ -194,16 +197,10 @@ def ocn_modelvsobs(config, field):

if buildObsClimatologies:
dsObs = xr.open_mfdataset(obsFileName)
dsTimeSlice = dsObs.sel(time=slice(timeStart, timeEnd))

# The following line converts from DASK to numpy to supress an odd
# warning that doesn't influence the figure output
dsTimeSlice.SSS.values

monthlyClimatology = dsTimeSlice.groupby('time.month').mean('time')

# Rename the observation data for code compactness
dsObs = monthlyClimatology.transpose('month', 'lat', 'lon')
dsObs.rename({'time': 'Time'}, inplace=True)
dsObs = dsObs.transpose('Time', 'lat', 'lon')
dsObs = dsObs.sel(Time=slice(timeStart, timeEnd))
dsObs['month'] = dsObs['Time.month']

obsFieldName = 'SSS'

Expand All @@ -225,8 +222,6 @@ def ocn_modelvsobs(config, field):
changed, startYear, endYear = \
climatology.update_start_end_year(ds, config, calendar)

monthlyClimatology = climatology.compute_monthly_climatology(ds, calendar)

mpasMappingFileName = climatology.write_mpas_mapping_file(
config=config, meshFileName=restartFileName)

Expand All @@ -253,8 +248,8 @@ def ocn_modelvsobs(config, field):
monthNames=months)

if overwriteMpasClimatology or not os.path.exists(climatologyFileName):
seasonalClimatology = climatology.compute_seasonal_climatology(
monthlyClimatology, monthValues, field)
seasonalClimatology = climatology.compute_climatology(
ds, monthValues, calendar)
# write out the climatology so we can interpolate it with
# interpolate.remap
seasonalClimatology.to_netcdf(climatologyFileName)
Expand Down Expand Up @@ -283,7 +278,7 @@ def ocn_modelvsobs(config, field):
if (overwriteObsClimatology or
(not os.path.exists(climatologyFileName) and
not os.path.exists(regriddedFileName))):
seasonalClimatology = climatology.compute_seasonal_climatology(
seasonalClimatology = climatology.compute_climatology(
dsObs, monthValues, obsFieldName)
# Either we want to overwite files or neither the climatology
# nor its regridded counterpart exist. Write out the
Expand Down
45 changes: 23 additions & 22 deletions mpas_analysis/sea_ice/modelvsobs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
Last Modified
-------------
03/23/2017
03/29/2017
"""

import os
Expand All @@ -34,8 +34,7 @@
from ..shared.generalized_reader.generalized_reader \
import open_multifile_dataset

from ..shared.timekeeping.utility import get_simulation_start_time, \
days_to_datetime
from ..shared.timekeeping.utility import get_simulation_start_time

from .utility import setup_sea_ice_task

Expand Down Expand Up @@ -65,7 +64,7 @@ def seaice_modelvsobs(config, streamMap=None, variableMap=None):
Last Modified
-------------
03/23/2017
03/29/2017
"""

# perform common setup for the task
Expand Down Expand Up @@ -128,20 +127,15 @@ def seaice_modelvsobs(config, streamMap=None, variableMap=None):
changed, startYear, endYear = \
climatology.update_start_end_year(ds, config, calendar)

monthlyClimatology = climatology.compute_monthly_climatology(ds, calendar)

mpasMappingFileName = climatology.write_mpas_mapping_file(
config=config, meshFileName=restartFileName)

_compute_and_plot_concentration(config, monthlyClimatology,
mpasMappingFileName)
_compute_and_plot_concentration(config, ds, mpasMappingFileName, calendar)

_compute_and_plot_thickness(config, monthlyClimatology,
mpasMappingFileName)
_compute_and_plot_thickness(config, ds, mpasMappingFileName, calendar)


def _compute_and_plot_concentration(config, monthlyClimatology,
mpasMappingFileName):
def _compute_and_plot_concentration(config, ds, mpasMappingFileName, calendar):
"""
Given a config file, monthly climatology on the mpas grid, and the data
necessary to perform horizontal interpolation to a comparison grid,
Expand All @@ -152,18 +146,22 @@ def _compute_and_plot_concentration(config, monthlyClimatology,
----------
config : an instance of MpasConfigParser
monthlyClimatology : an xarray data set containing a monthly climatology
ds : ``xarray.Dataset`` object
an xarray data set from which to compute climatologies
mpasMappingFileName : The name of a mapping file used to perform
interpolation of MPAS model results
calendar: ``{'gregorian', 'gregorian_noleap'}``
The name of one of the calendars supported by MPAS cores
Authors
-------
Xylar Asay-Davis, Milena Veneziani
Last Modified
-------------
03/03/2017
03/29/2017
"""

print " Make ice concentration plots..."
Expand Down Expand Up @@ -230,8 +228,8 @@ def _compute_and_plot_concentration(config, monthlyClimatology,
monthNames=months)

if overwriteMpasClimatology or not os.path.exists(climatologyFileName):
seasonalClimatology = climatology.compute_seasonal_climatology(
monthlyClimatology, monthValues, field)
seasonalClimatology = climatology.compute_climatology(
ds, monthValues, calendar)
# write out the climatology so we can interpolate it with
# interpolate.remap
seasonalClimatology.to_netcdf(climatologyFileName)
Expand Down Expand Up @@ -329,8 +327,7 @@ def _compute_and_plot_concentration(config, monthlyClimatology,
cbarlabel="fraction")


def _compute_and_plot_thickness(config, monthlyClimatology,
mpasMappingFileName):
def _compute_and_plot_thickness(config, ds, mpasMappingFileName, calendar):
"""
Given a config file, monthly climatology on the mpas grid, and the data
necessary to perform horizontal interpolation to a comparison grid,
Expand All @@ -341,18 +338,22 @@ def _compute_and_plot_thickness(config, monthlyClimatology,
----------
config : an instance of MpasConfigParser
monthlyClimatology : an xarray data set containing a monthly climatology
ds : ``xarray.Dataset`` object
an xarray data set from which to compute climatologies
mpasMappingFileName : The name of a mapping file used to perform
interpolation of MPAS model results
calendar: ``{'gregorian', 'gregorian_noleap'}``
The name of one of the calendars supported by MPAS cores
Authors
-------
Xylar Asay-Davis, Milena Veneziani
Last Modified
-------------
03/03/2017
03/29/2017
"""

print " Make ice thickness plots..."
Expand Down Expand Up @@ -411,8 +412,8 @@ def _compute_and_plot_thickness(config, monthlyClimatology,
monthNames=months)

if overwriteMpasClimatology or not os.path.exists(climatologyFileName):
seasonalClimatology = climatology.compute_seasonal_climatology(
monthlyClimatology, monthValues, field)
seasonalClimatology = climatology.compute_climatology(
ds, monthValues, calendar)
# write out the climatology so we can interpolate it with
# interpolate.remap. Set _FillValue so ncremap doesn't produce
# an error
Expand Down
Loading

0 comments on commit fefb7aa

Please sign in to comment.