Commit 2b56e664 authored by Lee de Mora's avatar Lee de Mora

Initial commit of all the files used to make the plots in gmd-2015-135

parents
class C2Chl:
"""class for relationships of carbon to chlorophyll following
Sathyendranath et al. 2009.
The class includes the loglog functions of the from log10(C)=m+p*log10(Chl) and their transformation C=m*Chl^p as well as the carbon to chlorophyll ratio C2Chl=m*Chl^(p-1)."""
def __init__(self,m,p):
self.loglog= lambda x:self._loglog(x,m,p)
self.C= lambda x:self._C(x,m,p)
self.C2Chl= lambda x:self._C2Chl(x,m,p)
_loglog = lambda self,lchl,m,p: m+p*lchl
_C = lambda self,chl,m,p: 10**m*chl**p
_C2Chl = lambda self,chl,m,p: 10**m*chl**(p-1)
def __call__(self,lchl):
return self.loglog(lchl)
#Sathyendranath HPLC (C_P):
SH = C2Chl(1.9,.65)
#Sathyendranath Turner (C_P):
ST = C2Chl(1.81,.63)
#Buck Turner:
B = C2Chl(1.92,.69)
#Sathyendranath HPLC (C_T):
SH_T = C2Chl(2.25,0.48)
#Sathyendranath Turner (C_T):
ST_T = C2Chl(2.19,0.45)
#
# Copyright 2015, Plymouth Marine Laboratory
#
# This file is part of the GMD-2015-135 toolkit, associated with the Geoscientific
# Model Development paper: gmd-2015-135,
# Submitted on 26 Jun 2015
# The role of ecosystem function and emergent relationships in the assessment of global marine ecosystem models: a case study with ERSEM
# L. de Mora, M. Butenschön, and J. I. Allen.
# Geosci. Model Dev. Discuss., 8, 6095-6141, 2015
# www.geosci-model-dev-discuss.net/8/6095/2015/
# doi:10.5194/gmdd-8-6095-2015
#
# The GMD-2015-135 toolkit is free software: you can redistribute it and/or modify it
# under the terms of the Revised Berkeley Software Distribution (BSD) 3-clause license.
# The GMD-2015-135 toolkit is distributed in the hope that it will be useful, but
# without any warranty; without even the implied warranty of merchantability
# or fitness for a particular purpose. See the revised BSD license for more details.
# You should have received a copy of the revised BSD license along with ukesm-validation.
# If not, see <http://opensource.org/licenses/BSD-3-Clause>.
#
# Address:
# Plymouth Marine Laboratory
# Prospect Place, The Hoe
# Plymouth, PL1 3DH, UK
#
# Email:
# ledm@pml.ac.uk
#
This toolkit was used to produce the plots from the paper:
The role of ecosystem function and emergent relationships in the assessment of global marine ecosystem models: a case study with ERSEM
L. de Mora, M. Butenschön, and J. I. Allen.
Geosci. Model Dev. Discuss., 8, 6095-6141, 2015
www.geosci-model-dev-discuss.net/8/6095/2015/
doi:10.5194/gmdd-8-6095-2015
Please note the following caveats:
1. This toolkit does not include model data.
2. This toolkit does not include any pre-processing tools.
3. This toolkit is supplied as is, and is may not be updated with patches or upgrades.
4. These tools are tailored to the iMarNet simulation of the ERSEM model,
in the ORCA1 global model grid, accessing data on the PML local network,
so changes to the code WILL be required to apply them in any other environment.
These tools require the followin python dictionaries:
matplotlib
netCDF4
numpy
scipy
shelve
pandas
They also require the following gitlab tools:
communitystructure
C2Chl
This diff is collapsed.
This diff is collapsed.
""" Toolkit for calculating the communicty structure.
It should work for single values and for array.
"""
from numpy import zeros,logspace,log10,exp,min,max,arange,clip
from comstrucFit import comstrucFit
# Taka's CSparameterisations:
pcclip = lambda a: clip(a*100.,0.,100.)
microHirata2011 = lambda c: clip(100./(0.9117 + exp(-2.733 * log10(c) + 0.4003)), 0., 100. )
picoHirata2011 = lambda c: clip(100.*(1./(-1.*(0.1529 + exp(1.0306 * log10(c) -1.5576))) - 1.8587*log10(c)+2.9954), 0., 100.)
nanoHirata2011 = lambda c: clip(100. - microHirata2011(c) - picoHirata2011(c), 0., 100.)
piconanoHirata2011 = lambda c: clip(100. - microHirata2011(c), 0., 100.)
CSparameterisations = ['hirata2011',"devred2011","brewin2012","brewin2011a","brewin2010", "brewin2014", "brewin2015","brewin2015b"]
#note:
# "brewin2015" is the corrected version
# "brewin2015b" is the original version that included a bug.
def continuumModel( chl,pft, fit="brewin2010",):
if fit=="brewin2010":
Cmpn = 1.057
Spn = 0.851
Cmp = 0.107
Sp = 6.801
CmpmSpm = 0.9
CmpSp = 0.728
if fit=="brewin2011a":
Cmpn = 0.775
Spn = 1.152
Cmp = 0.146
Sp = 5.118
CmpmSpm = 0.8933
CmpSp = 0.747
if fit=="brewin2012":
Cmpn = 0.937
Spn = 1.033
Cmp = 0.170
Sp = 4.804
CmpmSpm = 0.968
CmpSp = 0.817
if fit=="devred2011":
Cmpn = 0.546
Spn = 1.830
Cmp = 0.148
Sp = 6.765
CmpmSpm = 1.
CmpSp = 1.
if fit=="brewin2014":
Cmpn = 2.611
Spn = 0.364
Cmp = 0.730
Sp = 1.047
CmpmSpm = 0.951
CmpSp = 0.763
if fit=="brewin2015b":
Cmpn = 0.79
Spn = 1.16
Cmp = 0.15
Sp = 5.0
CmpmSpm = 0.87
CmpSp = 0.75
if fit=="brewin2015":
Cmpn = 0.77
Spn = 1.22
Cmp = 0.13
Sp = 6.15
CmpmSpm = 0.9394
CmpSp = 0.7995
p = Cmp *(1. - exp(-Sp *chl))
if pft=='pico': return p
pn = Cmpn*(1. - exp(-Spn*chl))
if pft=='piconano': return pn
if pft=='nano': return pn - p
if pft=='micro': return chl - pn
return False
def chlPercent(chl, pft,fit='Hirata 2011',):
""" Takes Chl in [mg Chl / m^3], a pft, and a fit, then returns the %.
"""
p= pft.lower()
fit = fit.lower().replace(' ','')
pfts = ['pico','nano','micro','piconano']
if p not in pfts:
print "chlPercent:\tERROR:\tWrong functional name", p, 'not in ', pfts
return False
if fit not in CSparameterisations:
print "chlPercent:\tERROR:\tWrong fit name", fit, 'not in ', CSparameterisations
return False
if fit == 'hirata2011':
if p == 'micro': return microHirata2011(chl)
if p == 'pico': return picoHirata2011(chl)
if p == 'nano': return nanoHirata2011(chl)
if p == 'piconano': return nanoHirata2011(chl)+picoHirata2011(chl)
return False
else:
return clip(100.*continuumModel( chl,pft, fit=fit)/chl,0.,100.)
def test():
""" Test the term chlPercent
"""
from itertoolsmodule import product
pfts = ['pico','nano','micro','piconano']
#CSparameterisations = ['Hirata2011',"Devred2011","Brewin2012","Brewin2011a","Brewin2010", "Brewin 2014",]
chls = [0.01,0.1,0.5,1.,1.5,2.,10.]
for f in CSparameterisations:
for p,c in product(pfts,chls):
print f,p,c,':\t',chlPercent(c, p,f)
from scipy.optimize import curve_fit
from numpy import ma,exp, arange, logspace,isnan, isinf, sqrt, diag
class comstrucFit:
""" This code produces a fit to data, using least squares regression.
Typically run in the following way:
fy = comstrucFit(totalchl,planktonFunctionTypePercent, pft='pico').plotForX(chlAxis)
"""
def Cp(self, tchl, Cmp, Sp ): return ma.array(Cmp * (1. - exp(-Sp *tchl))/tchl) #pico
def Cpn(self,tchl, Cmpn, Spn): return ma.array(Cmpn *(1. - exp(-Spn *tchl))/tchl) #piconano
#def pnmp(self,tchl, Cmp, Sp, Cmpn, Spn): return self.pn(tchl, Cmpn, Spn) - self.p(tchl, Cmp, Sp)#nano
#def cmpn(self,tchl, Cmpn, Spn): return tchl - self.pn( tchl, Cmpn, Spn) #micro
def microHirata2011(self,tchl, a0,a1,a2): return ma.aray(1./(a0 + exp(a1 * log10(tchl) + a2)))
def picoHirata2011( self,tchl, a0,a1,a2,a3,a4): return ma.array(1./(-1.*(a0 + exp(a1 * log10(c) + a2))) + a3*log10(c)+ a4)
# def microHirata2011 = lambda c: 1./(0.912 + exp(-2.733 * log10(c) + 0.4)))
# def picoHirata2011 = lambda c: (1./(-1.*(0.153 + exp(1.031 * log10(c) -1.558))) - 1.860*log10(c)+2.995))
# def nanoHirata2011 = lambda c: 1. - microHirata2011(c) - picoHirata2011(c))
#def piconanoHirata2011 = lambda c: 1. - microHirata2011(c))
def __init__(self,tchl, pftfrac, pft='pico',chlRange=[0,0],pcRange=[-4.,104.],fitTypes ='Continuum',):
if chlRange[0]==chlRange[1]:
chlRange = [tchl.min(),tchl.max()]
print "Community Structure Fit:\tSetting range:", chlRange
self.chlRange =chlRange
self.pcRange = pcRange
self.fitTypes = fitTypes
self.pft = pft
tchl = ma.array(tchl)
pftfrac = ma.array(pftfrac)
tchl = ma.masked_outside(tchl, chlRange[0], chlRange[1])
pftfrac = ma.masked_outside(pftfrac, pcRange[0], pcRange[1])
tchl = ma.masked_where(tchl.mask + pftfrac.mask, tchl )
pftfrac = ma.masked_where(tchl.mask + pftfrac.mask, pftfrac)
self.tchl = tchl.compressed()
self.pftfrac = pftfrac.compressed()
if fitTypes == 'runningMean':
self.tchl, self.pftfrac = self.doRunningMean()
if len(self.tchl)!=len(self.pftfrac):
print "Community Structure Fit:\tERROR:\tTotal Chl and ", self.pft," % do not match", len(self.tchl),'!=',len(self.pftfrac)
assert False
if self.fitTypes == 'Continuum':
if self.pft=='pico': self.func = self.Cp
if self.pft=='piconano':self.func = self.Cpn
if self.fitTypes == 'Hirata':
print "WARNING: This is not tested."
if self.pft=='pico': self.func = self.picoHirata2011
#if self.pft=='piconano':self.func = self.piconanoHirata2011
#if self.pft=='nano': self.func = self.nanoHirata2011
if self.pft=='micro': self.func = self.microHirata2011
#if self.pft not in ['pico', 'piconano',]:#'micro', 'nano',]:
# print "Community Structure Fit:\t not performing fit to ",self.pft
if 'NoFit' not in [fitTypes,]: self.fit()
def doRunningMean(self,nbins=1001):
""" Bin the two 1D arrays, to simplyfy the fit.
"""
print "Running Mean:", self.chlRange[0], ma.log10(self.chlRange[0]), self.chlRange[1], ma.log10(self.chlRange[1])
x_axis = logspace(ma.log10(self.chlRange[0]),ma.log10(self.chlRange[1]), nbins)
xbins = zip(x_axis[:-1],x_axis[1:])
data = {}
for x in xbins:
y = ma.masked_where((self.tchl<=x[0])+(self.tchl>x[1]), self.pftfrac).compressed()
x = ma.mean(x)
if y is ma.masked or not len(y):
data[x] = ma.masked
continue
#data[x] = ma.masked
print x, y.mean(), len(y)
if len(y)<len(self.tchl)/1000.:
data[x] = ma.masked
else: data[x] = ma.mean(y)
tchl = ma.array([ i for i in sorted(data.keys())])
pftfrac = ma.array([data[i] for i in sorted(data.keys())])
return tchl,pftfrac
#print "tchl:", self.tchl[0], self.tchl.mean(), self.tchl.min(), self.tchl.max()
#print "pftfrac:", self.pftfrac[0], self.pftfrac.mean(), self.pftfrac.min(), self.pftfrac.max()
def fit(self):
print "Community Structure Fit:\t performing fit",self.pft, self.fitTypes
# try: popt, pcov = curve_fit(self.func, self.tchl, self.pftfrac,)
try:
popt, pcov , infodict, errmsg, ier= curve_fit(self.func, self.tchl, self.pftfrac,full_output = True)
except:
print "Unable to solve and find a fit"
popt = ma.masked_all(2)
pcov = ma.masked_all((2,2))
infodict = []
errmsg = []
ier = []
print "Fitting: popt:", popt
print "Fitting: pcov:", pcov
print "Fitting: infodict:", infodict
print "Fitting: errmsg:", errmsg
print "Fitting: ier:", ier
#if self.pft=='pico': popt, pcov = curve_fit(self.p, self.tchl, self.pftfrac,)
#if self.pft=='piconano':popt, pcov = curve_fit(self.pn, self.tchl, self.pftfrac,)
#if self.pft=='nano': popt, pcov = curve_fit(self.pnmp, self.tchl, self.pftfrac,)
#if self.pft=='micro': popt, pcov = curve_fit(self.cmpn, self.tchl, self.pftfrac,)
print "Community Structure Fit:\t fit results:",popt, pcov
if ma.is_masked(pcov) or isnan(pcov).any() or isinf(pcov).any():
print "Community Structure Fit:\tWARNING:\tfit results contains a masked Value:",popt, pcov
self.popt , self.pcov = popt, pcov
if self.pft=='pico':
self.Cmp = self.popt[0]
self.Sp = self.popt[1]
self.Cmpn = ma.masked
self.Spn = ma.masked
if self.pft=='piconano':
self.Cmp = ma.masked
self.Sp = ma.masked
self.Cmpn = self.popt[0]
self.Spn = self.popt[1]
self.perr = sqrt(diag(pcov))
return popt, pcov
def plotForX(self,x):
if self.pft not in ['pico', 'piconano','micro', 'nano',]:
print "Community Structure Fit:\t not performing fit to ",self.pft
z = ma.zeros(len(x))
return ma.masked_where(z==0,z)
if ma.is_masked(self.pcov) or isnan(self.pcov).any() or isinf(self.pcov).any():
print "Community Structure Fit:\tWARNING:\tfit results contains a masked Value:",self.popt, self.pcov
z = ma.zeros(len(x))
return ma.masked_where(z==0,z)
print "Community Structure Fit:\t plot for x",self.pft,self.popt
y=[]
x = ma.array(x)
if x.min() < self.chlRange[0] or x.max() > self.chlRange[1]:
print "Community Structure Fit:\tWarning:\tplot for x ", self.pft," x ouside Chl Range", [x.min(),x,max()], '<',self.chlRange
#if:
# print "Community Structure Fit:\tWarning:\tplot for x ", self.pft," x > Chl Range", x.max(), '>',self.chlRange
if self.pft in ['pico', 'piconano','micro']: y = self.func(x, self.popt[0],self.popt[1])
if self.pft in ['nano',]: y = self.func(x, self.popt[0],self.popt[1],self.popt[2],self.popt[3])
#if self.pft=='pico':
#if self.pft=='piconano':y = self.pn(x, self.popt[0],self.popt[1])
#if self.pft=='nano': y = self.p(x, self.popt[0],self.popt[1])
#if self.pft=='micro': y = self.p(x, self.popt[0],self.popt[1])
x = ma.array(x)
y = ma.array(y)
print x, y
if len(y) != len(x):
print "Community Structure Fit:\tERROR:\tplot for x ", self.pft," % do not match", len(y),'!=',len(x)
assert False
if not len(y):
print "Community Structure Fit:\tERROR:\tCould not produce y in plotForX"
else: return y
#if self.pft=='piconano': popt, pcov = curve_fit(self.ppn, self.tchl, self.pftfrac,)
#if self.pft=='nano': popt, pcov = curve_fit(self.ppnmp, self.tchl, self.pftfrac,)
assert False
if __name__=="__main__":
#test()
#def test():
from ncdfView import ncdfView
fn = "/data/euryale7/scratch/ledm/iMarNet/xhonp/MEANS/xhonpo_20011231m01P.nc"
print "loading",fn
nc = ncdfView(fn, Quiet=True)
chl = nc('Chl3')[0,0]
tchl= nc('Chl1')[0,0]+nc('Chl2')[0,0]+nc('Chl3')[0,0]+nc('Chl4')[0,0]
pftfrac = chl/tchl
pft = 'pico'
popt, pcov = comstrucFit(tchl, pftfrac, pft=pft)
print 'The end.'
This diff is collapsed.
from os.path import exists
from os import mkdir, makedirs
from ncdfView import ncdfView
from glob import glob
from string import join
from mpl_toolkits.basemap import maskoceans
from netCDF4 import Dataset
from shelve import open as shOpen
from subprocess import Popen
from datetime import date,datetime, timedelta
from scipy.stats import linregress
from ICESmask import ICESmask
from matplotlib import pyplot
from pickle import dump, load
from numpy import meshgrid, fromfile, array, concatenate, where,delete
from numpy.ma import masked_where, zeros, arange, array as marray, masked
from DataPaths import DataPaths
from numpyXtns import point_inside_polygon
#--------------------------------------------------
# folder, lastWord
#--------------------------------------------------
def folder(name):
if type(name) == type(['a','b','c']):
name=join(name,'/')
if name[-1] != '/':
name = name+'/'
if exists(name) is False:
makedirs(name)
print 'makedirs ', name
return name
def lastWord(a, separator='/',debug=False):
#returns the final word of a string.
while a[-1] == separator: a = a[:-1]
ncount, count = 0,1
while ncount != count:
ncount = count
count = a[count:].find('/')+1+count
#if ncount == count : break
if debug: print 'final word:', a[count:]
return a[count:]
def addStraightLineFit(ax, x,y,showtext=True, addOneToOne=False,extent = [0,0,0,0]):
b1, b0, rValue, pValue, stdErr = getLinRegText(ax, x, y, showtext =showtext)
if extent == [0,0,0,0]:
fx = arange(x.min(), x.max(), (x.max()-x.min())/20.)
fy =[b0 + b1*a for a in fx]
else:
minv = min(extent)
maxv = max(extent)
fx = array(arange(minv, maxv, (maxv-minv)/1000.))
fy = array([b0 + b1*a for a in fx])
fx = masked_where((fx<minv) + (fy < minv) + (fx>maxv) + (fy > maxv), fx)
fy = masked_where((fx<minv) + (fy < minv) + (fx>maxv) + (fy > maxv), fy)
pyplot.plot(fx,fy, 'k')
if addOneToOne: pyplot.plot(fx,fx, 'k--')
#xstep = (x.max()-x.min())/40.
#ystep = (y.max()-y.min())/40.
#pyplot.axis([x.min()-xstep, x.max()+xstep, y.min()-ystep, y.max()+ystep])
def getLinRegText(ax, x, y, showtext=True):
x = [a for a in x if (a is masked)==False]
y = [a for a in y if (a is masked)==False]
beta1, beta0, rValue, pValue, stdErr = linregress(x, y)
thetext = r'$\^\beta_0$ = '+strRound(beta0) \
+ '\n'+r'$\^\beta_1$ = '+strRound(beta1) \
+ '\nR = '+ strRound(rValue) \
+ '\nP = '+strRound(pValue) \
+ '\nN = '+str(int(len(x)))
#+ '\n'+r'$\epsilon$ = ' + strRound(stdErr) \
if showtext: pyplot.text(0.04, 0.96,thetext ,
horizontalalignment='left',
verticalalignment='top',
transform = ax.transAxes)
return beta1, beta0, rValue, pValue, stdErr
def getIndex(x, arr):
"""getIndex: takes a depth and an array of box boundaries, returns depth index
if -1, coordinates are outside grid."""
x = float(x)
if ((arr[:,0]-x)*(arr[:,1]-x)).min() >0:return -1
return ((arr[:,0]-x)*(arr[:,1]-x)).argmin()
class AutoVivification(dict):
"""Implementation of perl's autovivification feature."""
def __getitem__(self, item):
try: return dict.__getitem__(self, item)
except KeyError:
value = self[item] = type(self)()
return value
def mnStr(month):
mn = '%02d' % month
return mn
def getLogTicks(mi,ma,quiet=False):
# this takes a min and max and produces a sensible log scale to use with logarithmic colorbars.
# cbar.set_ticks(n)
# cbar.set_ticklabels(n)
if mi*ma<0.:
print "Log scale should not be below Zero and above Zero."
assert False
n = [pow(10.,n) for n in xrange(-28,28) ]
while n[0] < mi: n = delete(n,0)
while n[-1] > ma: n = delete(n,-1)
if not quiet: print n, mi, ma
return n
This diff is collapsed.
#!/usr/bin/python
# -*- coding: latin-1 -*-
from sys import argv
from os import listdir,getcwd
from numpy import zeros
from netCDF4 import Dataset
from netCDF4 import num2date as n2d
from numpyXtns import spread
try:
from matplotlib.colors import LogNorm
from mpl_toolkits.basemap import Basemap
except:
print "No plotting capability loaded..."
try:
from xtraPlots import pmDarkMap,pmLightMap,pmMap,hovmoeller,worldPlot,discreteColors,chlMap,yearMap,colorWheel,zooMap,SMap,TMap,spmMap,pdfMap,mapMerc,mapOrtho,mapPlot,worldCylPlot,worldRegPlot,Faces,chlMapFun
except:
print "No extra plotting objects loaded..."
class ncdfView(Dataset):
def __enter__(self,):
return self
def __init__(self,filename,Mask=True,Quiet=False):
if not Quiet:
print 'Opening nCDF-file %s ...' % filename
Dataset.__init__(self,filename,'r')
if Mask:
for var in self.variables.values():
var.set_auto_maskandscale(True)
if not Quiet:
print self
def lon(self):
Lon=None
for key,var in self.variables.items():
try:
if var.units=='degrees_east':
print "...extracting longitude..."
try:
lon=self.variables[key+'_bnds'][:]
if lon.shape[1]==2:
Lon=zeros(lon.shape[0]+1)
Lon[:-1]=lon[:,0]
Lon[-1]=lon[-1,1]
print "bnd"
else:
Lon=lon
except:
Lon=var[:]
except:
pass
return Lon
def lat(self):
Lat=None
for key,var in self.variables.items():
try:
if var.units=='degrees_north':
print "...extracting latitude..."
try:
lat=self.variables[key+'_bnds'][:]
if lat.shape[1]==2:
Lat=zeros(lat.shape[0]+1)
Lat[:-1]=lat[:,0]
Lat[-1]=lat[-1,1]
print "bnd"
else:
Lat=lat
except:
Lat=var[:]
except:
pass
return Lat
def time(self):
Time=None
for key,Var in self.variables.items():
if key in ('time','days','hours','minutes','seconds'):
print "...found ",key,"..."
Time=Var[:]
else:
try:
Date=n2d(0,Var.units)
Time=Var[:]
print "...found ",key,"..."
except:
pass
return Time
def dates(self):
Dates=None
for key,Var in self.variables.items():
if key=='time':
print "...found ",key,"..."
n2d(0,Var.units)
try:
Dates=n2d(0,Var.units)
print "...found ",key,Var.units,"..."
print "Reference date of time variable: ",Dates
Dates=n2d(Var[:],Var.units)
except:
pass
return Dates
def __call__(self,varStr,Squeeze=False,Object=False):
if Object:
return self.variables[varStr]
else:
if Squeeze:
a=self.variables[varStr][:].squeeze()