Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change climatologies to use xtime_start and xtime_end #165

Merged
merged 4 commits into from
Apr 3, 2017
Merged

Change climatologies to use xtime_start and xtime_end #165

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
30e96d6
Add startTime and endTime coords in mpas_xarray
xylar Mar 29, 2017
8663690
Update climatologies to use start and end times
xylar Mar 29, 2017
c8b0d13
Add CI testing for climatology
xylar Mar 31, 2017
d65c74e
White space
xylar Mar 31, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion mpas_analysis/ocean/meridional_overturning_circulation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,8 @@
from ..shared.analysis_task import setup_task
from ..shared.climatology import climatology

from ..shared.constants import constants


def moc_streamfunction(config): # {{{
"""
Expand Down Expand Up @@ -304,7 +306,10 @@ def _compute_moc_climo_postprocess(config, runStreams, variableMap, calendar,
dictClimo['endYearClimo'] = endYear

# Compute annual climatology
annualClimatology = ds.mean('Time')
annualClimatology = \
climatology.compute_climatology(ds,
constants.monthDictionary['ANN'],
calendar)

# Convert to numpy arrays
# (can result in a memory error for large array size)
Expand Down
25 changes: 15 additions & 10 deletions mpas_analysis/ocean/nino34_index.py
View file
Open in desktop
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,9 @@ 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)
regionSST = ds.avgSurfaceTemperature.isel(nOceanRegions=regionIndex)
nino34 = compute_nino34_index(regionSST, calendar)

print ' Computing NINO3.4 power spectra...'
f, spectra, conf99, conf95, redNoise = compute_nino34_spectra(nino34)
Expand Down Expand Up @@ -122,7 +123,7 @@ def compute_nino34_index(regionSST, calendar): # {{{

Parameters
----------
regionSST : xarray.dataArray object
regionSST : xarray.DataArray object
values of SST in the nino region

calendar: {'gregorian', 'gregorian_noleap'}
Expand All @@ -138,18 +139,22 @@ def compute_nino34_index(regionSST, calendar): # {{{

Last Modified
-------------
03/23/2017
03/30/2017
"""

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


# add 'month' data array so we can group by month below.
regionSST = climatology.add_months_and_days_in_month(regionSST, 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(regionSST)

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

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


def compute_nino34_spectra(nino34Index): # {{{
Expand Down
55 changes: 25 additions & 30 deletions mpas_analysis/ocean/ocean_modelvsobs.py
View file
Open in desktop
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.coords['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')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@xylar, why is a transpose used here?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@pwolfram, the transpose is because the observational data is not in the same order as MPAS output.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@pwolfram, what's your concern about the transpose?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Transpose are notorious for being a bad idea, e.g., they typically indicate things weren't set up correctly and they are expensive to execute. It doesn't sound like it is a huge problem here but I just wanted to double-check.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are pretty small data sets and we've been transposing them the whole time. I just revised the code a bit.

dsObs = dsObs.sel(Time=slice(timeStart, timeEnd))
dsObs.coords['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.coords['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,8 +278,8 @@ def ocn_modelvsobs(config, field):
if (overwriteObsClimatology or
(not os.path.exists(climatologyFileName) and
not os.path.exists(regriddedFileName))):
seasonalClimatology = climatology.compute_seasonal_climatology(
dsObs, monthValues, obsFieldName)
seasonalClimatology = climatology.compute_climatology(
dsObs, monthValues)
# Either we want to overwite files or neither the climatology
# nor its regridded counterpart exist. Write out the
# climatology so we can interpolate it with interpolate.remap
Expand Down
45 changes: 23 additions & 22 deletions mpas_analysis/sea_ice/modelvsobs.py
View file
Open in desktop
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