diff --git a/gridded/da/analysis/siteseries.py b/gridded/da/analysis/siteseries.py
index c1169b15f9d2bd58d72cfad6cae10cc3ec14e463..f86f92024f2ea2c2f52a0f3fbbc444e6632ef873 100755
--- a/gridded/da/analysis/siteseries.py
+++ b/gridded/da/analysis/siteseries.py
@@ -11,6 +11,9 @@ File created on 23 Dec 2008.
import sys
sys.path.append('../../')
+import matplotlib
+matplotlib.use('pdf')
+
import sys
import os
from da.tools.general import create_dirs
@@ -22,6 +25,7 @@ import numpy as np
import datetime as dt
import da.tools.io4 as io
import logging
+import copy
from da.analysis.summarize_obs import nice_lon, nice_lat, nice_alt
import Image
import urllib2
@@ -47,7 +51,7 @@ fontsize = 9
figures that were created are stamped and saved properly
"""
-def site_timeseries(dacycle):
+def site_timeseries(analysisdir):
"""***************************************************************************************
Call example:
@@ -58,13 +62,14 @@ def site_timeseries(dacycle):
#
# Create directories if needed
#
- mrdir = os.path.join(dacycle['dir.analysis'], 'timeseries_molefractions')
+
+ mrdir = os.path.join(analysisdir, 'timeseries_molefractions')
if not os.path.exists(mrdir):
create_dirs(mrdir)
#
# Make a dictionary of available sites and the NetCDF file names associated with them
#
- sitelist = os.listdir(os.path.join(dacycle['dir.analysis'], 'data_molefractions'))
+ sitelist = os.listdir(os.path.join(analysisdir, 'data_molefractions'))
sitelist = [f for f in sitelist if f.endswith('.nc')]
#
# Loop over site and sitefiles
@@ -75,64 +80,68 @@ def site_timeseries(dacycle):
#
# Create filename and extract site codes for the sites that had day/night data separated
#
- filename = os.path.join(dacycle['dir.analysis'], 'data_molefractions', sitefile)
- #
- # Create a figure instance to hold plot
- #
- fig = plt.figure(1, figsize=(15, 6,))#,frameon=False)
- #
- # Make plot
- #
- fig = timevssite_new(fig, filename)
- #
- # Save image
- #
+ filename = os.path.join(analysisdir, 'data_molefractions', sitefile)
saveas = os.path.join(mrdir, sitefile[:-3] + '_timeseries')
- fig.savefig(saveas, dpi=100)
- #
- plt.close(fig)
+
+ if not os.path.exists(saveas+'.pdf'):
+ #
+ # Create a figure instance to hold plot
+ #
+ fig = plt.figure(1, figsize=(15, 6,))#,frameon=False)
+ #
+ # Make plot
+ #
+ fig = timevssite_new(fig, filename)
+ #
+ # Save image
+ #
+ fig.savefig(saveas, dpi=100)
+ #
+ plt.close(fig)
#
# Residuals next
#
logging.debug('Making residuals figures for %s: ' % sitefile)
- #
- # Create a figure instance to hold plot
- #
- fig = plt.figure(1, figsize=(15, 6,))#,frameon=False)
- #
- # Make plot
- #
- fig = residuals_new(fig, filename)
- #
- # Save image
- #
saveas = os.path.join(mrdir, sitefile[:-3] + '_residuals')
- fig.savefig(saveas, dpi=100)
- #
- # next in loop over sites
- #
- plt.close(fig)
+ if not os.path.exists(saveas+'.pdf'):
+ #
+ # Create a figure instance to hold plot
+ #
+ fig = plt.figure(1, figsize=(15, 6,))#,frameon=False)
+ #
+ # Make plot
+ #
+ fig = residuals_new(fig, filename)
+ #
+ # Save image
+ #
+ fig.savefig(saveas, dpi=100)
+ #
+ # next in loop over sites
+ #
+ plt.close(fig)
#
# histograms next
#
logging.debug('Making histograms figures for %s: ' % sitefile)
- #
- # Create a figure instance to hold plot
- #
- fig = plt.figure(1, figsize=(15, 6,))#,frameon=False)
- #
- # Make plot
- #
- fig = timehistograms_new(fig, filename)
- #
- # Save image
- #
saveas = os.path.join(mrdir, sitefile[:-3] + '_histograms')
- fig.savefig(saveas, dpi=100)
- #
- # next in loop over sites
- #
- plt.close(fig)
+ if not os.path.exists(saveas+'.pdf'):
+ #
+ # Create a figure instance to hold plot
+ #
+ fig = plt.figure(1, figsize=(15, 6,))#,frameon=False)
+ #
+ # Make plot
+ #
+ fig = timehistograms_new(fig, filename)
+ #
+ # Save image
+ #
+ fig.savefig(saveas, dpi=100)
+ #
+ # next in loop over sites
+ #
+ plt.close(fig)
def timehistograms_new(fig, infile):
"""
@@ -146,11 +155,21 @@ def timehistograms_new(fig, infile):
# Get data
#
f = io.CT_CDF(infile, 'read')
+ species = f.dataset_parameter
+ if species == 'co2':
+ molefac=1e6
+ units = '$\mu$mol mol$^{-1}$'
+ species = "CO$_2$"
+ if species == 'co2c13':
+ molefac=1.0
+ units = 'permil'
+ species = "$\delta^{13}$C"
+
date = f.get_variable('time')
- obs = f.get_variable('value') * 1e6
- mdm = f.get_variable('modeldatamismatch') * 1e6
- hphtr = f.get_variable('totalmolefractionvariance_forecast') * 1e6 * 1e6
- simulated = f.get_variable('modelsamplesmean_forecast') * 1e6
+ obs = f.get_variable('value') * molefac
+ mdm = f.get_variable('modeldatamismatch') * molefac
+ hphtr = f.get_variable('totalmolefractionvariance_forecast') * molefac * molefac
+ simulated = f.get_variable('modelsamplesmean_forecast') * molefac
flags = f.get_variable('flag_forecast')
longsitestring = f.site_name + ', ' + f.site_country
@@ -278,16 +297,14 @@ def timehistograms_new(fig, infile):
ax.xaxis.set_ticks_position('bottom')
- ax.set_xlabel('[$\mu$mol mol$^{-1}$]')
+ ax.set_xlabel('[%s]'%units)
#
# All custom titles and auxiliary info are placed onto the figure directly (fig.text) in relative coordinates
#
- fig.text(0.5, 0.02, 'Simulated - Observed CO$_2$ ($\mu$mol/mol)\nData from %s to %s' %
- (pydates[0].strftime('%d-%b-%Y'), pydates[-1].strftime('%d-%b-%Y')),
- horizontalalignment='center', fontsize=fontsize - 1)
+ fig.text(0.5, 0.02, 'Simulated - Observed %s [%s]\nData from %s to %s' %(species,units,pydates[0].strftime('%d-%b-%Y'), pydates[-1].strftime('%d-%b-%Y')), horizontalalignment='center', fontsize=fontsize - 1)
#fig.text(0.75,0.02,'Simulated - Observed\n CO$_2$ ($\mu$mol/mol)',horizontalalignment='center',fontsize=fontsize)
- fig.text(0.5, 0.35, 'model-data\nmismatch:\n%.2f $\mu$mol/mol' % sel_mdm[0], horizontalalignment='center', fontsize=fontsize, color='green')
+ fig.text(0.5, 0.35, 'model-data\nmismatch:\n%.2f %s' % (sel_mdm[0], units), horizontalalignment='center', fontsize=fontsize, color='green')
fig.text(0.12, 0.75, 'NH Summer\n(Jun-Sep)', horizontalalignment='center', fontsize=fontsize)
fig.text(0.62, 0.75, 'NH Winter\n(Nov-Apr)', horizontalalignment='center', fontsize=fontsize)
#
@@ -334,10 +351,19 @@ def timevssite_new(fig, infile):
# Get data
#
f = io.CT_CDF(infile, 'read')
+ species = f.dataset_parameter
+ if species == 'co2':
+ molefac=1e6
+ units = '$\mu$mol mol$^{-1}$'
+ species = "CO$_2$"
+ if species == 'co2c13':
+ molefac=1.0
+ units = 'permil'
+ species = "$\delta^{13}$C"
date = f.get_variable('time')
- obs = f.get_variable('value') * 1e6
- mdm = f.get_variable('modeldatamismatch') * 1e6
- simulated = f.get_variable('modelsamplesmean') * 1e6
+ obs = f.get_variable('value') * molefac
+ mdm = f.get_variable('modeldatamismatch') * molefac
+ simulated = f.get_variable('modelsamplesmean') * molefac
flags = f.get_variable('flag_forecast')
longsitestring = f.site_name + ', ' + f.site_country
@@ -417,7 +443,7 @@ def timevssite_new(fig, infile):
#
# Location and format of xticks
#
- ax1.xaxis.set_minor_locator(pltdt.MonthLocator())
+ ax1.xaxis.set_minor_locator(pltdt.MonthLocator([1,7],bymonthday=1))
ax1.xaxis.set_minor_formatter(pltdt.DateFormatter('%b'))
ax1.xaxis.set_major_locator(pltdt.YearLocator())
ax1.xaxis.set_major_formatter(pltdt.DateFormatter('\n\n%Y'))
@@ -445,7 +471,7 @@ def timevssite_new(fig, infile):
#xtitle='Time'
#ax1.set_xlabel(xtitle, fontsize=fontsize) # label x axis
- ax1.set_ylabel(r"CO$_2$ [ppm]", fontsize=fontsize) # label y-axis
+ ax1.set_ylabel(r"%s [%s]"% (species,units), fontsize=fontsize) # label y-axis
#
# Axes 2
@@ -537,11 +563,20 @@ def residuals_new(fig, infile):
# Get data
#
f = io.CT_CDF(infile, 'read')
+ species = f.dataset_parameter
+ if species == 'co2':
+ molefac=1e6
+ units = '$\mu$mol mol$^{-1}$'
+ species = "CO$_2$"
+ if species == 'co2c13':
+ molefac=1.0
+ units = 'permil'
+ species = "$\delta^{13}$C"
date = f.get_variable('time')
- obs = f.get_variable('value') * 1e6
- mdm = f.get_variable('modeldatamismatch') * 1e6
- hphtr = f.get_variable('totalmolefractionvariance_forecast') * 1e6 * 1e6
- simulated = f.get_variable('modelsamplesmean_forecast') * 1e6
+ obs = f.get_variable('value') * molefac
+ mdm = f.get_variable('modeldatamismatch') * molefac
+ simulated = f.get_variable('modelsamplesmean_forecast') * molefac
+ hphtr = f.get_variable('totalmolefractionvariance_forecast') * molefac * molefac
flags = f.get_variable('flag_forecast')
longsitestring = f.site_name + ', ' + f.site_country
@@ -672,7 +707,7 @@ def residuals_new(fig, infile):
#
# Location and format of xticks
#
- ax1.xaxis.set_minor_locator(pltdt.MonthLocator())
+ ax1.xaxis.set_minor_locator(pltdt.MonthLocator([1,7],bymonthday=1))
ax1.xaxis.set_minor_formatter(pltdt.DateFormatter('%b'))
ax1.xaxis.set_major_locator(pltdt.YearLocator())
ax1.xaxis.set_major_formatter(pltdt.DateFormatter('\n\n%Y'))
@@ -699,7 +734,7 @@ def residuals_new(fig, infile):
#xtitle='Time'
#ax1.set_xlabel(xtitle, fontsize=fontsize) # label x axis
- ax1.set_ylabel(r"CO$_2$ [ppm]", fontsize=fontsize) # label y-axis
+ ax1.set_ylabel(r"%s [%s]"%(species,units), fontsize=fontsize) # label y-axis
#
# Title
#
@@ -743,27 +778,13 @@ def residuals_new(fig, infile):
if __name__ == '__main__': # started as script
- from da.tools.initexit import CycleControl
- from da.carbondioxide.dasystem import CO2DaSystem
- from da.carbondioxide.statevector import CO2StateVector
-
sys.path.append('../../')
logging.root.setLevel(logging.DEBUG)
- dacycle = CycleControl(args={'rc':'../../ctdas-od-gfed2-glb6x4-obspack-full.rc'})
- dacycle.initialize()
- dacycle.parse_times()
-
- dasystem = CO2DaSystem('../rc/carbontracker_ct09_opf.rc')
- dasystem.initialize()
-
- dacycle.dasystem = dasystem
+ analysisdir = "/Storage/CO2/peters/CTDAS-OD-GFED2-GLB6x4-ObsPack-full-karolina/analysis/"
- q = site_timeseries(dacycle)
- #q=site_residuals(dacycle)
- #q=site_statistics(rundat)
- #q=site_histograms(rundat)
+ site_timeseries(analysisdir)
sys.exit(0)