analysis_level3_dms.py

#!/usr/bin/ipython 

#
# Copyright 2015, Plymouth Marine Laboratory
#
# This file is part of the bgc-val library.
#
# bgc-val is free software: you can redistribute it and/or modify it
# under the terms of the Revised Berkeley Software Distribution (BSD) 3-clause license. 

# bgc-val 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 bgc-val.
# 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
#
"""
.. module:: analysis_level3_dms
   :platform: Unix
   :synopsis: A script to produce a level 3 analysis for DMS.

.. moduleauthor:: Lee de Mora <ledm@pml.ac.uk>

"""


#####	
# Load Standard Python modules:
from sys import argv,exit
from os.path import exists
from calendar import month_name
from socket import gethostname
from netCDF4 import Dataset
from glob import glob
from scipy.interpolate import interp1d
import numpy as np
import os,sys
from getpass import getuser

#####	
# Load specific local code:
import UKESMpython as ukp
from timeseries import timeseriesAnalysis
from timeseries import profileAnalysis
from timeseries import timeseriesPlots as tsp 

#####
# User defined set of paths pointing towards the datasets.
import paths	


medusaCoords 	= {'t':'time_counter', 'z':'deptht', 'lat': 'nav_lat',  'lon': 'nav_lon',   'cal': '360_day',}	# model doesn't need time dict.
dmsCoords	= {'t':'time',    'z':'depth',  'lat':'Latitude',  'lon': 'Longitude','cal': 'standard','tdict':ukp.tdicts['ZeroToZero']}	



def analysis_dms(jobID=''):
	annual = True
	
	analysisDict = {}
	imagedir	 = ukp.folder(paths.imagedir+'/'+jobID+'/Level3/DMS')
	shelvedir 	= ukp.folder(paths.shelvedir+'/'+jobID+'/Level3/DMS')

	regionList	= ['Global', 'ignoreInlandSeas',
	  		'SouthernOcean','ArcticOcean',
			'Equator10', 'Remainder',
			'NorthernSubpolarAtlantic','NorthernSubpolarPacific',
			]
  	metricList = ['mean',]			
	dataD = {}		
	modeldataD = {}


	def listModelDataFiles(jobID, filekey, datafolder, annual):
		if annual:
			return sorted(glob(datafolder+jobID+"/"+jobID+"o_1y_*_"+filekey+".nc"))
		else:
			return sorted(glob(datafolder+jobID+"/"+jobID+"o_1m_*_"+filekey+".nc"))
					        			



	#####
	# A time series analysis for the DMS fields.
	
	for name in ['DMS_ARAN','DMS_ANDR','DMS_SIMO','DMS_HALL']:


		dmsfiles = listModelDataFiles(jobID, 'diad_T', paths.ModelFolder_pref, annual)
		if name == 'DMS_ARAN':	
			analysisDict['modelFiles'] = dmsfiles
		else:
			analysisDict['modelFiles'] = ukp.listFiles(dmsfiles, want=100, listType='backloaded',first=30,last=10)
		
		if annual:		
			analysisDict['dataFile'] = paths.DMSDir+'DMSclim_mean.nc'
		else:	analysisDict['dataFile'] = ''
			
		analysisDict['modelcoords'] 	= medusaCoords 	
		analysisDict['datacoords'] 	= dmsCoords
	
		analysisDict['modeldetails'] 	= {'name': name, 'vars':['DMS_ARAN',], 'convert': ukp.mul1000000,'units':'umol/m3'}
		analysisDict['datadetails']  	= {'name': name, 'vars':['DMS',],      'convert': ukp.NoChange,'units':'umol/m3'}
	
		analysisDict['layers'] 		= ['layerless',]
		analysisDict['regions'] 	= regionList	
		analysisDict['metrics']		= metricList

		analysisDict['datasource'] 	= 'Lana'
		analysisDict['model']		= 'MEDUSA'

		analysisDict['modelgrid']	= 'eORCA1'
		analysisDict['gridFile']	= paths.orcaGridfn
		analysisDict['Dimensions']	= 2	
		
		
		
		tsa = timeseriesAnalysis(
		        analysisDict['modelFiles'],
		        analysisDict['dataFile'],
		        dataType        = name,
		        modelcoords     = analysisDict['modelcoords'],
		        modeldetails    = analysisDict['modeldetails'],
		        datacoords      = analysisDict['datacoords'],
		        datadetails     = analysisDict['datadetails'],
		        datasource      = analysisDict['datasource'],
		        model           = analysisDict['model'],
		        jobID           = jobID,
		        layers          = analysisDict['layers'],
		        regions         = analysisDict['regions'],
		        metrics         = analysisDict['metrics'],
		        workingDir      = shelvedir,
		        imageDir        = imagedir,
		        grid            = analysisDict['modelgrid'],
		        gridFile        = analysisDict['gridFile'],
		        clean           = False,
		)

		dataD[name] = tsa.dataD
		modeldataD[name] = tsa.modeldataD


	#####
	# Prepare a time series comparison of the four DMS types.
	
	timesD  = {}
	arrD	= {}
	
	region = 'Global'
	for name in dataD.keys():
		try:mdata = modeldataD[name][(region, 'layerless', 'mean')]
		except: continue	
		
		timesD[name] 	= sorted(mdata.keys())
		arrD[name]	= [mdata[t] for t in timesD[name]]
		

	colours = {'DMS_ARAN':'red', 'DMS_SIMO':'orange','DMS_ANDR':'blue','DMS_HALL':'purple', }
	title = 'DMS ' + region		
	
	for ls in ['Both','smoothed','movingaverage']:
		tsp.multitimeseries(
			timesD, 		# model times (in floats)
			arrD,			# model time series
			data 	= -999,		# in situ data distribution
			title 	= title,
			filename=ukp.folder(imagedir)+'DMS_'+region+'_'+ls+'.png',
			units = '',
			plotStyle 	= 'Together',
			lineStyle	= ls,
			colours		= colours,
		)
			
			

def main():
	try:	jobID = argv[1]
	except:	
		jobID = "u-ab749"

	if 'debug' in argv[1:]:	suite = 'debug'
		
	
	analysis_dms(jobID =jobID, )#clean=1)			
	#if suite == 'all':
	#analysis_timeseries(jobID =jobID,analysisSuite='FullDepth', z_component = 'FullDepth',)#clean=1)  

if __name__=="__main__":
	main()