From 32be3ea2b45a8a498e8f38c2ae6b51ffa9e64061 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 7 Dec 2016 17:47:39 +0100 Subject: [PATCH] Add variable and stream maps to sea-ice analysis Also, make sea-ice scripts PEP8 compliant. --- mpas_analysis/sea_ice/modelvsobs.py | 273 +++++++++++-------- mpas_analysis/sea_ice/timeseries.py | 211 +++++++------- mpas_analysis/sea_ice/variable_stream_map.py | 16 ++ run_analysis.py | 10 +- 4 files changed, 298 insertions(+), 212 deletions(-) create mode 100644 mpas_analysis/sea_ice/variable_stream_map.py diff --git a/mpas_analysis/sea_ice/modelvsobs.py b/mpas_analysis/sea_ice/modelvsobs.py index 6f576d038..bbd2d0332 100644 --- a/mpas_analysis/sea_ice/modelvsobs.py +++ b/mpas_analysis/sea_ice/modelvsobs.py @@ -8,22 +8,23 @@ import numpy.ma as ma import xarray as xr import datetime -#import calendar + from netCDF4 import Dataset as netcdf_dataset -from ..shared.mpas_xarray.mpas_xarray import preprocess_mpas, remove_repeated_time_index +from ..shared.mpas_xarray.mpas_xarray import preprocess_mpas, \ + remove_repeated_time_index from ..shared.plot.plotting import plot_polar_comparison from ..shared.io import StreamsFile -from ..shared.io.utility import paths -def seaice_modelvsobs(config): + +def seaice_modelvsobs(config, streamMap=None, variableMap=None): """ Performs analysis of sea-ice properties by comparing with previous model results and/or observations. Author: Xylar Asay-Davis, Milena Veneziani - Last Modified: 10/27/2016 + Last Modified: 12/07/2016 """ # read parameters from config file @@ -36,99 +37,127 @@ def seaice_modelvsobs(config): # reading only those that are between the start and end dates startDate = config.get('time', 'climo_start_date') endDate = config.get('time', 'climo_end_date') - infiles = streams.readpath('timeSeriesStatsMonthlyOutput', - startDate=startDate, endDate=endDate) - print 'Reading files {} through {}'.format(infiles[0],infiles[-1]) + streamName = streams.find_stream(streamMap['timeSeriesStats']) + infiles = streams.readpath(streamName, startDate=startDate, + endDate=endDate) + print 'Reading files {} through {}'.format(infiles[0], infiles[-1]) - plots_dir = config.get('paths','plots_dir') - obsdir = config.get('paths','obs_seaicedir') + plots_dir = config.get('paths', 'plots_dir') + obsdir = config.get('paths', 'obs_seaicedir') - casename = config.get('case','casename') + casename = config.get('case', 'casename') - remapfile = config.get('data','mpas_remapfile') - climodir = config.get('data','mpas_climodir') + remapfile = config.get('data', 'mpas_remapfile') + climodir = config.get('data', 'mpas_climodir') - climo_yr1 = config.getint('time','climo_yr1') - climo_yr2 = config.getint('time','climo_yr2') - yr_offset = config.getint('time','yr_offset') + climo_yr1 = config.getint('time', 'climo_yr1') + climo_yr2 = config.getint('time', 'climo_yr2') + yr_offset = config.getint('time', 'yr_offset') - #climodir = "%s/%s" % (climodir,casename) - climodir_regridded = "%s/mpas_regridded" % climodir - if not os.path.isdir("%s" % climodir): + # climodir = "{}/{}".format(climodir, casename) + climodir_regridded = "{}/mpas_regridded".format(climodir) + if not os.path.isdir(climodir): print "\nClimatology directory does not exist. Create it...\n" - os.mkdir("%s" % climodir) - if not os.path.isdir("%s" % climodir_regridded): + os.mkdir(climodir) + if not os.path.isdir(climodir_regridded): print "\nRegridded directory does not exist. Create it...\n" - os.mkdir("%s" % climodir_regridded) + os.mkdir(climodir_regridded) print indir print climodir # Model climo (output) filenames climofiles = {} - climofiles['winNH'] = "mpas-cice_climo.years%04d-%04d.jfm.nc" % (climo_yr1,climo_yr2) - climofiles['sumNH'] = "mpas-cice_climo.years%04d-%04d.jas.nc" % (climo_yr1,climo_yr2) - climofiles['winSH'] = "mpas-cice_climo.years%04d-%04d.djf.nc" % (climo_yr1,climo_yr2) - climofiles['sumSH'] = "mpas-cice_climo.years%04d-%04d.jja.nc" % (climo_yr1,climo_yr2) - climofiles['on'] = "mpas-cice_climo.years%04d-%04d.on.nc" % (climo_yr1,climo_yr2) - climofiles['fm'] = "mpas-cice_climo.years%04d-%04d.fm.nc" % (climo_yr1,climo_yr2) + climofiles['winNH'] = "mpas-cice_climo.years{:04d}-{:04d}.jfm.nc".format( + climo_yr1, climo_yr2) + climofiles['sumNH'] = "mpas-cice_climo.years{:04d}-{:04d}.jas.nc".format( + climo_yr1, climo_yr2) + climofiles['winSH'] = "mpas-cice_climo.years{:04d}-{:04d}.djf.nc".format( + climo_yr1, climo_yr2) + climofiles['sumSH'] = "mpas-cice_climo.years{:04d}-{:04d}.jja.nc".format( + climo_yr1, climo_yr2) + climofiles['on'] = "mpas-cice_climo.years{:04d}-{:04d}.on.nc".format( + climo_yr1, climo_yr2) + climofiles['fm'] = "mpas-cice_climo.years{:04d}-{:04d}.fm.nc".format( + climo_yr1, climo_yr2) # make a dictionary of the months in each climotology monthsInClim = {} - monthsInClim['winNH'] = [1,2,3] - monthsInClim['sumNH'] = [7,8,9] - monthsInClim['winSH'] = [12,1,2] - monthsInClim['sumSH'] = [6,7,8] - monthsInClim['on'] = [10,11] - monthsInClim['fm'] = [2,3] - + monthsInClim['winNH'] = [1, 2, 3] + monthsInClim['sumNH'] = [7, 8, 9] + monthsInClim['winSH'] = [12, 1, 2] + monthsInClim['sumSH'] = [6, 7, 8] + monthsInClim['on'] = [10, 11] + monthsInClim['fm'] = [2, 3] # Obs filenames obs_iceconc_filenames = {} - obs_iceconc_filenames['winNH_NASATeam'] = "%s/SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_NH_jfm.interp0.5x0.5.nc" % obsdir - obs_iceconc_filenames['sumNH_NASATeam'] = "%s/SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_NH_jas.interp0.5x0.5.nc" % obsdir - obs_iceconc_filenames['winSH_NASATeam'] = "%s/SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_SH_djf.interp0.5x0.5.nc" % obsdir - obs_iceconc_filenames['sumSH_NASATeam'] = "%s/SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_SH_jja.interp0.5x0.5.nc" % obsdir - obs_iceconc_filenames['winNH_Bootstrap'] = "%s/SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_NH_jfm.interp0.5x0.5.nc" % obsdir - obs_iceconc_filenames['sumNH_Bootstrap'] = "%s/SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_NH_jas.interp0.5x0.5.nc" % obsdir - obs_iceconc_filenames['winSH_Bootstrap'] = "%s/SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_SH_djf.interp0.5x0.5.nc" % obsdir - obs_iceconc_filenames['sumSH_Bootstrap'] = "%s/SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_SH_jja.interp0.5x0.5.nc" % obsdir + obs_iceconc_filenames['winNH_NASATeam'] = \ + "{}/SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_NH_" \ + "jfm.interp0.5x0.5.nc".format(obsdir) + obs_iceconc_filenames['sumNH_NASATeam'] = \ + "{}/SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_NH_" \ + "jas.interp0.5x0.5.nc".format(obsdir) + obs_iceconc_filenames['winSH_NASATeam'] = \ + "{}/SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_SH_" \ + "djf.interp0.5x0.5.nc".format(obsdir) + obs_iceconc_filenames['sumSH_NASATeam'] = \ + "{}/SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_SH_" \ + "jja.interp0.5x0.5.nc".format(obsdir) + obs_iceconc_filenames['winNH_Bootstrap'] = \ + "{}/SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_" \ + "NH_jfm.interp0.5x0.5.nc".format(obsdir) + obs_iceconc_filenames['sumNH_Bootstrap'] = \ + "{}/SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_" \ + "NH_jas.interp0.5x0.5.nc".format(obsdir) + obs_iceconc_filenames['winSH_Bootstrap'] = \ + "{}/SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_" \ + "SH_djf.interp0.5x0.5.nc".format(obsdir) + obs_iceconc_filenames['sumSH_Bootstrap'] = \ + "{}/SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_" \ + "SH_jja.interp0.5x0.5.nc".format(obsdir) obs_icethick_filenames = {} - obs_icethick_filenames['onNH'] = "%s/ICESat/ICESat_gridded_mean_thickness_NH_on.interp0.5x0.5.nc" % obsdir - obs_icethick_filenames['fmNH'] = "%s/ICESat/ICESat_gridded_mean_thickness_NH_fm.interp0.5x0.5.nc" % obsdir - obs_icethick_filenames['onSH'] = "%s/ICESat/ICESat_gridded_mean_thickness_SH_on.interp0.5x0.5.nc" % obsdir - obs_icethick_filenames['fmSH'] = "%s/ICESat/ICESat_gridded_mean_thickness_SH_fm.interp0.5x0.5.nc" % obsdir + obs_icethick_filenames['onNH'] = "{}/ICESat/ICESat_gridded_mean_" \ + "thickness_NH_on.interp0.5x0.5.nc".format(obsdir) + obs_icethick_filenames['fmNH'] = "{}/ICESat/ICESat_gridded_mean_" \ + "thickness_NH_fm.interp0.5x0.5.nc".format(obsdir) + obs_icethick_filenames['onSH'] = "{}/ICESat/ICESat_gridded_mean_" \ + "thickness_SH_on.interp0.5x0.5.nc".format(obsdir) + obs_icethick_filenames['fmSH'] = "{}/ICESat/ICESat_gridded_mean_" \ + "thickness_SH_fm.interp0.5x0.5.nc".format(obsdir) # Checks on directory/files existence: for climName in obs_iceconc_filenames: obs_filename = obs_iceconc_filenames[climName] - if os.path.isfile("%s" % obs_filename) != True: - raise SystemExit("Obs file %s not found. Exiting..." % obs_filename) + if not os.path.isfile(obs_filename): + raise SystemExit("Obs file {} not found. Exiting...".format( + obs_filename)) for climName in obs_icethick_filenames: obs_filename = obs_icethick_filenames[climName] - if os.path.isfile("%s" % obs_filename) != True: - raise SystemExit("Obs file %s not found. Exiting..." % obs_filename) - + if not os.path.isfile(obs_filename): + raise SystemExit("Obs file {} not found. Exiting...".format( + obs_filename)) # Load data print " Load sea-ice data..." - ds = xr.open_mfdataset(infiles, preprocess=lambda x: - preprocess_mpas(x, yearoffset=yr_offset, - timestr='timeSeriesStatsMonthly_avg_daysSinceStartOfSim_1', - onlyvars=['timeSeriesStatsMonthly_avg_iceAreaCell_1', - 'timeSeriesStatsMonthly_avg_iceVolumeCell_1'])) + ds = xr.open_mfdataset( + infiles, + preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset, + timestr='Time', + onlyvars=['iceAreaCell', + 'iceVolumeCell'], + variable_map=variableMap)) ds = remove_repeated_time_index(ds) # Compute climatologies (first motnhly and then seasonally) print " Compute seasonal climatologies..." - time_start = datetime.datetime(yr_offset+climo_yr1,1,1) - time_end = datetime.datetime(yr_offset+climo_yr2,12,31) - ds_tslice = ds.sel(Time=slice(time_start,time_end)) + time_start = datetime.datetime(yr_offset+climo_yr1, 1, 1) + time_end = datetime.datetime(yr_offset+climo_yr2, 12, 31) + ds_tslice = ds.sel(Time=slice(time_start, time_end)) # check that each year has 24 months (?) monthly_clim = ds_tslice.groupby('Time.month').mean('Time') - daysInMonth = [31,28,31,30,31,30,31,31,30,31,30,31] - monthLetters = ['J','F','M','A','M','J','J','A','S','O','N','D'] - + daysInMonth = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] + monthLetters = ['J', 'F', 'M', 'A', 'M', 'J', 'J', 'A', 'S', 'O', 'N', 'D'] clims = {} for climName in monthsInClim: @@ -148,9 +177,10 @@ def seaice_modelvsobs(config): print " Regrid fields to regular grid..." for climName in clims: # Save to netcdf files - outFileName = "%s/%s" % (climodir,climofiles[climName]) + outFileName = "{}/{}".format(climodir, climofiles[climName]) clims[climName].to_netcdf(outFileName) - args = ["ncremap", "-P", "mpas", "-i", outFileName, "-m", remapfile, "-O", climodir_regridded] + args = ["ncremap", "-P", "mpas", "-i", outFileName, "-m", remapfile, + "-O", climodir_regridded] try: subprocess.check_call(args) except subprocess.CalledProcessError, e: @@ -163,7 +193,7 @@ def seaice_modelvsobs(config): # interate over observations of sea-ice concentration first = True - for climName in ['winNH','winSH','sumNH','sumSH']: + for climName in ['winNH', 'winSH', 'sumNH', 'sumSH']: hemisphere = climName[-2:] season = climName[:-2] @@ -175,26 +205,29 @@ def seaice_modelvsobs(config): clevsModelObs = config.getExpression('seaice_modelvsobs', 'clevsModelObs_conc_{}'.format( season)) - cmap = plt.get_cmap(config.get('seaice_modelvsobs','cmapModelObs')) + cmap = plt.get_cmap(config.get('seaice_modelvsobs', 'cmapModelObs')) cmapIndices = config.getExpression('seaice_modelvsobs', 'cmapIndicesModelObs') - cmapModelObs = cols.ListedColormap(cmap(cmapIndices),"cmapModelObs") + cmapModelObs = cols.ListedColormap(cmap(cmapIndices), "cmapModelObs") clevsDiff = config.getExpression('seaice_modelvsobs', 'clevsDiff_conc_{}'.format(season)) - cmap = plt.get_cmap(config.get('seaice_modelvsobs','cmapDiff')) + cmap = plt.get_cmap(config.get('seaice_modelvsobs', 'cmapDiff')) cmapIndices = config.getExpression('seaice_modelvsobs', 'cmapIndicesDiff') - cmapDiff = cols.ListedColormap(cmap(cmapIndices),"cmapDiff") + cmapDiff = cols.ListedColormap(cmap(cmapIndices), "cmapDiff") - lon0 = config.getfloat('seaice_modelvsobs','lon0_%s'%hemisphere) - latmin = config.getfloat('seaice_modelvsobs','latmin_%s'%hemisphere) + lon0 = config.getfloat('seaice_modelvsobs', + 'lon0_{}'.format(hemisphere)) + latmin = config.getfloat('seaice_modelvsobs', + 'latmin_{}'.format(hemisphere)) # Load in sea-ice data # Model... # ice concentrations - f = netcdf_dataset("%s/%s" % (climodir_regridded,climofiles[climName]),mode='r') - iceconc = f.variables["timeSeriesStatsMonthly_avg_iceAreaCell_1"][:] + fileName = "{}/{}".format(climodir_regridded, climofiles[climName]) + f = netcdf_dataset(fileName, mode='r') + iceconc = f.variables["iceAreaCell"][:] if(first): lons = f.variables["lon"][:] lats = f.variables["lat"][:] @@ -206,9 +239,10 @@ def seaice_modelvsobs(config): # ...and observations # ice concentrations from NASATeam (or Bootstrap) algorithm - for obsName in ['NASATeam','Bootstrap']: + for obsName in ['NASATeam', 'Bootstrap']: - f = netcdf_dataset(obs_iceconc_filenames['%s_%s'%(climName, obsName)],mode='r') + fileName = obs_iceconc_filenames['{}_{}'.format(climName, obsName)] + f = netcdf_dataset(fileName, mode='r') obs_iceconc = f.variables["AICE"][:] f.close() @@ -219,6 +253,11 @@ def seaice_modelvsobs(config): monthsName.append(monthLetters[month-1]) monthsName = ''.join(monthsName) + title = "{} ({}, years {:04d}-{:04d})".format( + suptitle, monthsName, climo_yr1, climo_yr2) + fileout = "{}/iceconc{}{}_{}_{}_years{:04d}-{:04d}.png".format( + plots_dir, obsName, hemisphere, casename, monthsName, + climo_yr1, climo_yr2) plot_polar_comparison( config, Lons, @@ -230,29 +269,28 @@ def seaice_modelvsobs(config): clevsModelObs, cmapDiff, clevsDiff, - title = "%s (%s, years %04d-%04d)" % (suptitle, monthsName, climo_yr1, climo_yr2), - fileout = "%s/iceconc%s%s_%s_%s_years%04d-%04d.png" % (plots_dir, obsName, hemisphere, - casename, monthsName, climo_yr1, climo_yr2), - plotProjection = plotProjection, - latmin = latmin, - lon0 = lon0, - modelTitle = "%s" % casename, - obsTitle = "Observations (SSM/I %s)" % (obsName), - diffTitle = "Model-Observations", - cbarlabel = "%") - + title=title, + fileout=fileout, + plotProjection=plotProjection, + latmin=latmin, + lon0=lon0, + modelTitle=casename, + obsTitle="Observations (SSM/I {})".format(obsName), + diffTitle="Model-Observations", + cbarlabel="%") print " Make ice thickness plots..." # Plot Northern Hemisphere FM sea-ice thickness suptitle = "Ice thickness" # interate over observations of sea-ice thickness - for climName in ['fm','on']: + for climName in ['fm', 'on']: # Load in sea-ice data # Model... # ice concentrations - f = netcdf_dataset("%s/%s" % (climodir_regridded,climofiles[climName]),mode='r') - icethick = f.variables["timeSeriesStatsMonthly_avg_iceVolumeCell_1"][:] + fileName = "{}/{}".format(climodir_regridded, climofiles[climName]) + f = netcdf_dataset(fileName, mode='r') + icethick = f.variables["iceVolumeCell"][:] f.close() monthsName = [] @@ -260,46 +298,55 @@ def seaice_modelvsobs(config): monthsName.append(monthLetters[month-1]) monthsName = ''.join(monthsName) - for hemisphere in ['NH','SH']: + for hemisphere in ['NH', 'SH']: # ...and observations # ice concentrations from NASATeam (or Bootstrap) algorithm clevsModelObs = config.getExpression( 'seaice_modelvsobs', 'clevsModelObs_thick_{}'.format(hemisphere)) - cmap = plt.get_cmap(config.get('seaice_modelvsobs','cmapModelObs')) + cmap = plt.get_cmap(config.get('seaice_modelvsobs', + 'cmapModelObs')) cmapIndices = config.getExpression('seaice_modelvsobs', 'cmapIndicesModelObs') - cmapModelObs = cols.ListedColormap(cmap(cmapIndices),"cmapModelObs") + cmapModelObs = cols.ListedColormap(cmap(cmapIndices), + "cmapModelObs") clevsDiff = config.getExpression( 'seaice_modelvsobs', 'clevsDiff_thick_{}'.format(hemisphere)) - cmap = plt.get_cmap(config.get('seaice_modelvsobs','cmapDiff')) + cmap = plt.get_cmap(config.get('seaice_modelvsobs', 'cmapDiff')) cmapIndices = config.getExpression('seaice_modelvsobs', 'cmapIndicesDiff') - cmapDiff = cols.ListedColormap(cmap(cmapIndices),"cmapDiff") + cmapDiff = cols.ListedColormap(cmap(cmapIndices), "cmapDiff") - lon0 = config.getfloat('seaice_modelvsobs','lon0_%s'%hemisphere) - latmin = config.getfloat('seaice_modelvsobs','latmin_%s'%hemisphere) + lon0 = config.getfloat('seaice_modelvsobs', + 'lon0_{}'.format(hemisphere)) + latmin = config.getfloat('seaice_modelvsobs', + 'latmin_{}'.format(hemisphere)) - f = netcdf_dataset(obs_icethick_filenames['%s%s'%(climName,hemisphere)],mode='r') + fileName = obs_icethick_filenames['{}{}'.format(climName, + hemisphere)] + f = netcdf_dataset(fileName, mode='r') obs_icethick = f.variables["HI"][:] f.close() # Mask thickness fields - icethick[ icethick == 0 ] = ma.masked - obs_icethick = ma.masked_values(obs_icethick,0) + icethick[icethick == 0] = ma.masked + obs_icethick = ma.masked_values(obs_icethick, 0) if hemisphere == 'NH': # Obs thickness should be nan above 86 (ICESat data) - obs_icethick[ Lats > 86 ] = ma.masked + obs_icethick[Lats > 86] = ma.masked plotProjection = 'npstere' else: plotProjection = 'spstere' diff = icethick - obs_icethick - - + title = "{} ({}, years {:04d}-{:04d})".format(suptitle, monthsName, + climo_yr1, climo_yr2) + fileout = "{}/icethick{}_{}_{}_years{:04d}-{:04d}.png".format( + plots_dir, hemisphere, casename, monthsName, climo_yr1, + climo_yr2) plot_polar_comparison( config, Lons, @@ -311,14 +358,12 @@ def seaice_modelvsobs(config): clevsModelObs, cmapDiff, clevsDiff, - title = "%s (%s, years %04d-%04d)" % (suptitle, monthsName, climo_yr1, climo_yr2), - fileout = "%s/icethick%s_%s_%s_years%04d-%04d.png" % (plots_dir, hemisphere, - casename, monthsName, climo_yr1, climo_yr2), - plotProjection = plotProjection, - latmin = latmin, - lon0 = lon0, - modelTitle = "%s" % casename, - obsTitle = "Observations (ICESat)", - diffTitle = "Model-Observations", - cbarlabel = "m") - + title=title, + fileout=fileout, + plotProjection=plotProjection, + latmin=latmin, + lon0=lon0, + modelTitle=casename, + obsTitle="Observations (ICESat)", + diffTitle="Model-Observations", + cbarlabel="m") diff --git a/mpas_analysis/sea_ice/timeseries.py b/mpas_analysis/sea_ice/timeseries.py index 3ffac1168..ec9c20c1a 100644 --- a/mpas_analysis/sea_ice/timeseries.py +++ b/mpas_analysis/sea_ice/timeseries.py @@ -2,23 +2,23 @@ import xarray as xr import pandas as pd import datetime -import os.path -from ..shared.mpas_xarray.mpas_xarray import preprocess_mpas, remove_repeated_time_index +from ..shared.mpas_xarray.mpas_xarray import preprocess_mpas, \ + remove_repeated_time_index -from ..shared.plot.plotting import timeseries_analysis_plot +from ..shared.plot.plotting import timeseries_analysis_plot from ..shared.io import StreamsFile -from ..shared.io.utility import paths from ..shared.timekeeping.Date import Date -def seaice_timeseries(config): + +def seaice_timeseries(config, streamMap=None, variableMap=None): """ Performs analysis of time series of sea-ice properties. Author: Xylar Asay-Davis, Milena Veneziani - Last Modified: 11/28/2016 + Last Modified: 12/07/2016 """ # read parameters from config file @@ -31,54 +31,58 @@ def seaice_timeseries(config): # reading only those that are between the start and end dates startDate = config.get('time', 'timeseries_start_date') endDate = config.get('time', 'timeseries_end_date') - infiles = streams.readpath('timeSeriesStatsMonthlyOutput', - startDate=startDate, endDate=endDate) - print 'Reading files {} through {}'.format(infiles[0],infiles[-1]) - - varnames=['iceAreaCell','iceVolumeCell'] + streamName = streams.find_stream(streamMap['timeSeriesStats']) + infiles = streams.readpath(streamName, startDate=startDate, + endDate=endDate) + print 'Reading files {} through {}'.format(infiles[0], infiles[-1]) - plot_titles = {'iceAreaCell':'Sea-ice area', - 'iceVolumeCell':'Sea-ice volume', - 'iceThickness':'Sea-ice thickness'} + varnames = ['iceAreaCell', 'iceVolumeCell'] - units_dict = {'iceAreaCell':'[km$^2$]', - 'iceVolumeCell':'[10$^3$ km$^3$]', - 'iceThickness':'[m]'} + plot_titles = {'iceAreaCell': 'Sea-ice area', + 'iceVolumeCell': 'Sea-ice volume', + 'iceThickness': 'Sea-ice thickness'} - obs_filenames = {'iceAreaCell':[config.get('seaIceData','obs_iceareaNH'), - config.get('seaIceData','obs_iceareaSH')], - 'iceVolumeCell':[config.get('seaIceData','obs_icevolNH'), - config.get('seaIceData','obs_icevolSH')]} + units_dict = {'iceAreaCell': '[km$^2$]', + 'iceVolumeCell': '[10$^3$ km$^3$]', + 'iceThickness': '[m]'} + obs_filenames = { + 'iceAreaCell': [config.get('seaIceData', 'obs_iceareaNH'), + config.get('seaIceData', 'obs_iceareaSH')], + 'iceVolumeCell': [config.get('seaIceData', 'obs_icevolNH'), + config.get('seaIceData', 'obs_icevolSH')]} # Some plotting rules title_font_size = config.get('seaice_timeseries', 'title_font_size') - indir = config.get('paths','archive_dir_ocn') - meshfile = config.get('data','mpas_meshfile') + indir = config.get('paths', 'archive_dir_ocn') + meshfile = config.get('data', 'mpas_meshfile') - casename = config.get('case','casename') - ref_casename_v0 = config.get('case','ref_casename_v0') - indir_v0data = config.get('paths','ref_archive_v0_seaicedir') + casename = config.get('case', 'casename') + ref_casename_v0 = config.get('case', 'ref_casename_v0') + indir_v0data = config.get('paths', 'ref_archive_v0_seaicedir') - compare_with_obs = config.getboolean('seaice_timeseries','compare_with_obs') + compare_with_obs = config.getboolean('seaice_timeseries', + 'compare_with_obs') - plots_dir = config.get('paths','plots_dir') + plots_dir = config.get('paths', 'plots_dir') - yr_offset = config.getint('time','yr_offset') + yr_offset = config.getint('time', 'yr_offset') - N_movavg = config.getint('seaice_timeseries','N_movavg') + N_movavg = config.getint('seaice_timeseries', 'N_movavg') print " Load sea-ice data..." # Load mesh dsmesh = xr.open_dataset(meshfile) # Load data - ds = xr.open_mfdataset(infiles, preprocess=lambda x: - preprocess_mpas(x, yearoffset=yr_offset, - timestr='timeSeriesStatsMonthly_avg_daysSinceStartOfSim_1', - onlyvars=['timeSeriesStatsMonthly_avg_iceAreaCell_1', - 'timeSeriesStatsMonthly_avg_iceVolumeCell_1'])) + ds = xr.open_mfdataset( + infiles, + preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset, + timestr='Time', + onlyvars=['iceAreaCell', + 'iceVolumeCell'], + variable_map=variableMap)) ds = remove_repeated_time_index(ds) # convert the start and end dates to datetime objects using @@ -92,15 +96,17 @@ def seaice_timeseries(config): ds = ds.merge(dsmesh) year_start = (pd.to_datetime(ds.Time.min().values)).year - year_end = (pd.to_datetime(ds.Time.max().values)).year - time_start = datetime.datetime(year_start,1,1) - time_end = datetime.datetime(year_end,12,31) + year_end = (pd.to_datetime(ds.Time.max().values)).year + time_start = datetime.datetime(year_start, 1, 1) + time_end = datetime.datetime(year_end, 12, 31) if ref_casename_v0 != "None": - infiles_v0data = "".join([indir_v0data,'/icevol.',ref_casename_v0,'.year*.nc']) - ds_v0 = xr.open_mfdataset(infiles_v0data, preprocess=lambda x: - preprocess_mpas(x, yearoffset=yr_offset)) - ds_v0_tslice = ds_v0.sel(Time=slice(time_start,time_end)) + infiles_v0data = "{}/icevol.{}.year*.nc".format(indir_v0data, + ref_casename_v0) + ds_v0 = xr.open_mfdataset( + infiles_v0data, + preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) + ds_v0_tslice = ds_v0.sel(Time=slice(time_start, time_end)) # Make Northern and Southern Hemisphere partition: areaCell = ds.areaCell @@ -115,11 +121,11 @@ def seaice_timeseries(config): plot_title = plot_titles[varname] units = units_dict[varname] - print " Compute NH and SH time series of %s..."%(varname) + print " Compute NH and SH time series of {}...".format(varname) if varname == "iceThickCell": - varnamefull = "timeSeriesStatsMonthly_avg_iceVolumeCell_1" + varnamefull = "iceVolumeCell" else: - varnamefull = "".join(["timeSeriesStatsMonthly_avg_",varname,"_1"]) + varnamefull = varname var = ds[varnamefull] var_nh = var.where(ind_nh)*areaCell_nh @@ -132,15 +138,15 @@ def seaice_timeseries(config): if varname == "iceAreaCell": var_nh = var_nh.sum('nCells') var_sh = var_sh.sum('nCells') - var_nh = 1e-6*var_nh # m^2 to km^2 - var_sh = 1e-6*var_sh # m^2 to km^2 + var_nh = 1e-6*var_nh # m^2 to km^2 + var_sh = 1e-6*var_sh # m^2 to km^2 var_nh_iceext = 1e-6*var_nh_iceext.sum('nCells') var_sh_iceext = 1e-6*var_sh_iceext.sum('nCells') elif varname == "iceVolumeCell": var_nh = var_nh.sum('nCells') var_sh = var_sh.sum('nCells') - var_nh = 1e-3*1e-9*var_nh # m^3 to 10^3 km^3 - var_sh = 1e-3*1e-9*var_sh # m^3 to 10^3 km^3 + var_nh = 1e-3*1e-9*var_nh # m^3 to 10^3 km^3 + var_sh = 1e-3*1e-9*var_sh # m^3 to 10^3 km^3 else: var_nh = var_nh.mean('nCells')/areaCell_nh.mean('nCells') var_sh = var_sh.mean('nCells')/areaCell_sh.mean('nCells') @@ -150,67 +156,82 @@ def seaice_timeseries(config): xlabel = "Time [years]" if ref_casename_v0 != "None": - figname_nh = "%s/%sNH_%s_%s.png" % (plots_dir,varname,casename,ref_casename_v0) - figname_sh = "%s/%sSH_%s_%s.png" % (plots_dir,varname,casename,ref_casename_v0) + figname_nh = "{}/{}NH_{}_{}.png".format(plots_dir, varname, + casename, ref_casename_v0) + figname_sh = "{}/{}SH_{}_{}.png".format(plots_dir, varname, + casename, ref_casename_v0) else: - figname_nh = "%s/%sNH_%s.png" % (plots_dir,varname,casename) - figname_sh = "%s/%sSH_%s.png" % (plots_dir,varname,casename) + figname_nh = "{}/{}NH_{}.png".format(plots_dir, varname, casename) + figname_sh = "{}/{}SH_{}.png".format(plots_dir, varname, casename) - title_nh = "%s (NH), %s (r)" % (plot_title,casename) - title_sh = "%s (SH), %s (r)" % (plot_title,casename) + title_nh = "{} (NH), {} (r)".format(plot_title, casename) + title_sh = "{} (SH), {} (r)".format(plot_title, casename) if compare_with_obs: if varname == "iceAreaCell": - title_nh = "%s\nSSM/I observations, annual cycle (k)" % title_nh - title_sh = "%s\nSSM/I observations, annual cycle (k)" % title_sh + title_nh = \ + "{}\nSSM/I observations, annual cycle (k)".format(title_nh) + title_sh = \ + "{}\nSSM/I observations, annual cycle (k)".format(title_sh) elif varname == "iceVolumeCell": - title_nh = "%s\nPIOMAS, annual cycle (k)" % title_nh - title_sh = "%s\n" % title_sh + title_nh = "{}\nPIOMAS, annual cycle (k)".format(title_nh) + title_sh = "{}\n".format(title_sh) if ref_casename_v0 != "None": - title_nh = "%s\n %s (b)" % (title_nh,ref_casename_v0) - title_sh = "%s\n %s (b)" % (title_sh,ref_casename_v0) - + title_nh = "{}\n {} (b)".format(title_nh, ref_casename_v0) + title_sh = "{}\n {} (b)".format(title_sh, ref_casename_v0) if varname == "iceAreaCell": if compare_with_obs: - ds_obs = xr.open_mfdataset(obs_filenameNH, preprocess=lambda x: - preprocess_mpas(x, yearoffset=yr_offset)) + ds_obs = xr.open_mfdataset( + obs_filenameNH, + preprocess=lambda x: preprocess_mpas(x, + yearoffset=yr_offset)) ds_obs = remove_repeated_time_index(ds_obs) var_nh_obs = ds_obs.IceArea - var_nh_obs = replicate_cycle(var_nh,var_nh_obs) + var_nh_obs = replicate_cycle(var_nh, var_nh_obs) - ds_obs = xr.open_mfdataset(obs_filenameSH, preprocess=lambda x: - preprocess_mpas(x, yearoffset=yr_offset)) + ds_obs = xr.open_mfdataset( + obs_filenameSH, + preprocess=lambda x: preprocess_mpas(x, + yearoffset=yr_offset)) ds_obs = remove_repeated_time_index(ds_obs) var_sh_obs = ds_obs.IceArea - var_sh_obs = replicate_cycle(var_sh,var_sh_obs) + var_sh_obs = replicate_cycle(var_sh, var_sh_obs) if ref_casename_v0 != "None": - infiles_v0data = "".join([indir_v0data,'/icearea.',ref_casename_v0,'.year*.nc']) - ds_v0 = xr.open_mfdataset(infiles_v0data, preprocess=lambda x: - preprocess_mpas(x, yearoffset=yr_offset)) - ds_v0_tslice = ds_v0.sel(Time=slice(time_start,time_end)) + infiles_v0data = "{}/icearea.{}.year*.nc".format( + indir_v0data, ref_casename_v0) + ds_v0 = xr.open_mfdataset( + infiles_v0data, + preprocess=lambda x: preprocess_mpas(x, + yearoffset=yr_offset)) + ds_v0_tslice = ds_v0.sel(Time=slice(time_start, time_end)) var_nh_v0 = ds_v0_tslice.icearea_nh var_sh_v0 = ds_v0_tslice.icearea_sh elif varname == "iceVolumeCell": if compare_with_obs: - ds_obs = xr.open_mfdataset(obs_filenameNH, preprocess=lambda x: - preprocess_mpas(x, yearoffset=yr_offset)) + ds_obs = xr.open_mfdataset( + obs_filenameNH, + preprocess=lambda x: preprocess_mpas(x, + yearoffset=yr_offset)) ds_obs = remove_repeated_time_index(ds_obs) var_nh_obs = ds_obs.IceVol - var_nh_obs = replicate_cycle(var_nh,var_nh_obs) + var_nh_obs = replicate_cycle(var_nh, var_nh_obs) var_sh_obs = None if ref_casename_v0 != "None": - infiles_v0data = "".join([indir_v0data,'/icevol.',ref_casename_v0,'.year*.nc']) - ds_v0 = xr.open_mfdataset(infiles_v0data, preprocess=lambda x: - preprocess_mpas(x, yearoffset=yr_offset)) - ds_v0_tslice = ds_v0.sel(Time=slice(time_start,time_end)) + infiles_v0data = "{}/icevol.{}.year*.nc".format( + indir_v0data, ref_casename_v0) + ds_v0 = xr.open_mfdataset( + infiles_v0data, + preprocess=lambda x: preprocess_mpas(x, + yearoffset=yr_offset)) + ds_v0_tslice = ds_v0.sel(Time=slice(time_start, time_end)) var_nh_v0 = ds_v0_tslice.icevolume_nh var_sh_v0 = ds_v0_tslice.icevolume_sh @@ -248,39 +269,41 @@ def seaice_timeseries(config): title_font_size=title_font_size) else: # we will combine north and south onto a single graph - figname = "%s/%s.%s.png" % (plots_dir,casename,varname) - title = "%s, NH (r), SH (k)\n%s" % (plot_title,casename) + figname = "{}/{}.{}.png".format(plots_dir, casename, varname) + title = "{}, NH (r), SH (k)\n{}".format(plot_title, casename) timeseries_analysis_plot(config, [var_nh, var_sh], N_movavg, title, xlabel, units, figname, - lineStyles=['r-','k-'], + lineStyles=['r-', 'k-'], lineWidths=[1.2, 1.2], title_font_size=title_font_size) elif varname == "iceThickCell": - figname = "%s/%s.%s.png" % (plots_dir,casename,varname) - title = "%s NH (r), SH (k)\n%s" % (plot_title,casename) - timeseries_analysis_plot(config, [var_nh,var_sh], N_movavg, title, - xlabel, units, figname, - lineStyles=['r-','k-'], - lineWidths=[1.2, 1.2], - title_font_size=title_font_size) + figname = "{}/{}.{}.png".format(plots_dir, casename, varname) + title = "{} NH (r), SH (k)\n{}".format(plot_title, casename) + timeseries_analysis_plot(config, [var_nh, var_sh], N_movavg, title, + xlabel, units, figname, + lineStyles=['r-', 'k-'], + lineWidths=[1.2, 1.2], + title_font_size=title_font_size) else: - raise SystemExit("varname variable %s not supported for plotting" % varname) + raise ValueError( + "varname variable {} not supported for plotting".format( + varname)) -def replicate_cycle(ds,ds_toreplicate): +def replicate_cycle(ds, ds_toreplicate): dsshift = ds_toreplicate.copy() - shiftT = ((dsshift.Time.max() - dsshift.Time.min()) - + (dsshift.Time.isel(Time=1) - dsshift.Time.isel(Time=0))) + shiftT = ((dsshift.Time.max() - dsshift.Time.min()) + + (dsshift.Time.isel(Time=1) - dsshift.Time.isel(Time=0))) nT = np.ceil((ds.Time.max() - ds.Time.min())/shiftT) # replicate cycle: for i in np.arange(nT): dsnew = ds_toreplicate.copy() dsnew['Time'] = dsnew.Time + (i+1)*shiftT - dsshift = xr.concat([dsshift,dsnew], dim='Time') + dsshift = xr.concat([dsshift, dsnew], dim='Time') # constrict replicated ds_short to same time dimension as ds_long: dsshift = dsshift.sel(Time=ds.Time.values, method='nearest') return dsshift diff --git a/mpas_analysis/sea_ice/variable_stream_map.py b/mpas_analysis/sea_ice/variable_stream_map.py new file mode 100644 index 000000000..10fa93ab3 --- /dev/null +++ b/mpas_analysis/sea_ice/variable_stream_map.py @@ -0,0 +1,16 @@ +# mappings of stream names from various MPAS-SI versions to those in +# mpas_analysis +seaIceStreamMap = {'timeSeriesStats': ['timeSeriesStatsMonthlyOutput']} + + +seaIceVariableMap = {} +seaIceVariableMap['Time'] = \ + [['xtime_start', 'xtime_end'], + 'timeSeriesStatsMonthly_avg_daysSinceStartOfSim_1', + 'time_avg_daysSinceStartOfSim'] + +seaIceVariableMap['iceAreaCell'] = \ + ['timeSeriesStatsMonthly_avg_iceAreaCell_1'] + +seaIceVariableMap['iceVolumeCell'] = \ + ['timeSeriesStatsMonthly_avg_iceVolumeCell_1'] diff --git a/run_analysis.py b/run_analysis.py index b3ac15b53..864d77327 100644 --- a/run_analysis.py +++ b/run_analysis.py @@ -18,8 +18,8 @@ from mpas_analysis.ocean.variable_stream_map import oceanStreamMap, \ oceanVariableMap -# from mpas_analysis.sea_ice.variable_stream_map import seaIceStreamMap, \ -# seaIceVariableMap +from mpas_analysis.sea_ice.variable_stream_map import seaIceStreamMap, \ + seaIceVariableMap def path_existence(config, section, option, ignorestr=None): # {{{ @@ -142,14 +142,16 @@ def analysis(config): # {{{ print "" print "Plotting sea-ice area and volume time series..." from mpas_analysis.sea_ice.timeseries import seaice_timeseries - seaice_timeseries(config) + seaice_timeseries(config, streamMap=seaIceStreamMap, + variableMap=seaIceVariableMap) if config.getboolean('seaice_modelvsobs', 'generate'): print "" print "Plotting 2-d maps of sea-ice concentration and thickness " \ "climatologies..." from mpas_analysis.sea_ice.modelvsobs import seaice_modelvsobs - seaice_modelvsobs(config) + seaice_modelvsobs(config, streamMap=seaIceStreamMap, + variableMap=seaIceVariableMap) # GENERATE LAND-ICE DIAGNOSTICS