Commit 6dad6b69 authored by Peters, Wouter's avatar Peters, Wouter
Browse files

updated figure style

parent 11e7555a
......@@ -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)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment