fix previous codes and main work on the visualization and plot functi…

…ons (#11)

Modify wread.py.
Modify code pvalstats.py with plots functions.

postproc/comparevalww3board.py

@@ -52,8 +52,7 @@
 import matplotlib.colors as colors
 from matplotlib.colors import BoundaryNorm
 # import pickle
-import warnings
-warnings.filterwarnings("ignore")
+import warnings; warnings.filterwarnings("ignore")
 import wread
 palette = plt.cm.jet
 dpalette = plt.cm.RdBu_r

postproc/wread.py

@@ -45,7 +45,7 @@
 """
 
 import matplotlib
-# matplotlib.use('Agg')
+matplotlib.use('Agg')
 import time
 from time import strptime
 from calendar import timegm
@@ -225,11 +225,14 @@ def spec_ndbc(*args):
 	'''
 	Observations NDBC, wave spectrum, netcdf format
 	Input: file name (example: 46047w2016.nc)
-	Output: Values: time(datetime64),time(seconds since 1970),lat,lon; Arrays: freq,dfreq,pspec,dmspec,dpspec,r1spec,r2spec
+	Output: Values: time(datetime64),time(seconds since 1970),lat,lon; Arrays: freq,dfreq,pspec,dmspec,dpspec,dirspec
 	'''
-	if len(args) == 1:
+	deltatheta=np.int(10)
+	if len(args) >= 1:
 		fname=np.str(args[0])
-	elif len(args) > 1:
+	if len(args) >= 2:
+		deltatheta=np.int(args[1])
+	if len(args) > 2:
 		sys.exit(' Too many inputs')
 
 	try:
@@ -249,20 +252,27 @@ def spec_ndbc(*args):
 		r1spec = ds['wave_spectrum_r1'][:,:,0,0]
 		r2spec = ds['wave_spectrum_r2'][:,:,0,0]
 		ds.close(); del ds
-		# # DF in frequency (dfreq), https://www.ndbc.noaa.gov/wavespectra.shtml
+		# DF in frequency (dfreq), https://www.ndbc.noaa.gov/wavespectra.shtml
 		dfreq=np.zeros(47,'f')
 		dfreq[0]=0.010; dfreq[1:14]=0.005; dfreq[14:40]=0.010; dfreq[40::]=0.020
 		pspec=np.array(pspec*dfreq)
+		# Directional 2D Spectrum, https://www.ndbc.noaa.gov/measdes.shtml#swden , https://www.ndbc.noaa.gov/wavemeas.pdf
+		theta = np.array(np.arange(0,360+0.1,deltatheta))
+		# final directional wave spectrum (frequency X direction)
+		dirspec = np.zeros((btime.shape[0],freq.shape[0],theta.shape[0]),'f')
+		for t in range(0,btime.shape[0]):
+			dirspec[t,:,:] = np.array([pspec[t,:]]).T * (1/pi)*(0.5+  np.array([r1spec[t,:]]).T * cos(np.array( np.array([theta])-np.array([dmspec[t,:]]).T )*(pi/180)) 
+				+ np.array([r2spec[t,:]]).T*cos(2*np.array( np.array([theta]) - np.array([dpspec[t,:]]).T )*(pi/180)))
 
-	return bdate,btime,blat,blon,freq,dfreq,pspec,dmspec,dpspec,r1spec,r2spec
+	return bdate,btime,blat,blon,freq,dfreq,pspec,dmspec,dpspec,theta,dirspec
 
 
 # WAVEWATCH III spectra output, netcdf format
 def spec_ww3(*args):
 	'''
 	WAVEWATCH III, wave spectrum, netcdf format
 	Input: file name (example: ww3gefs.20160928_spec.nc), and station name (example: 41002)
-	Output: Values: time(datetime64),time(seconds since 1970),lat,lon; Arrays: freq,dfreq,pwst,d1sp,dire,dspec
+	Output: Values: time(datetime64),time(seconds since 1970),lat,lon; Arrays: freq,dfreq,pwst,d1sp,dire,dspec,wnds,wndd
 	'''
 	if len(args) == 2:
 		fname=np.str(args[0]); stname=np.str(args[1])

postproc/ww3fields.py

@@ -59,8 +59,7 @@
 from cartopy.util import add_cyclic_point
 from matplotlib import ticker
 # import pickle
-import warnings
-warnings.filterwarnings("ignore")
+import warnings; warnings.filterwarnings("ignore")
 
 palette = plt.cm.jet
 

postproc/ww3pointimeseries.py

@@ -43,8 +43,7 @@
 import matplotlib.pyplot as plt
 import sys
 import pandas as pd
-import warnings
-warnings.filterwarnings("ignore")
+import warnings; warnings.filterwarnings("ignore")
 
 # Input arguments
 if len(sys.argv) < 3 :

postproc/ww3pointspec.py

@@ -51,8 +51,7 @@
 from calendar import timegm
 import matplotlib.colors as colors
 import matplotlib.gridspec as gridspec
-import warnings
-warnings.filterwarnings("ignore")
+import warnings; warnings.filterwarnings("ignore")
 # Palette and colors for plotting the figures
 palette = plt.cm.gist_stern_r
 

validation/gridSatGlobal_Altimeter.py

@@ -52,8 +52,7 @@
 import time
 from calendar import timegm
 import sys
-import warnings
-warnings.filterwarnings("ignore")
+import warnings; warnings.filterwarnings("ignore")
 # netcdf format
 fnetcdf="NETCDF4"
 

validation/modelSat_collocation.py

@@ -27,7 +27,6 @@
    files/years in that directory will be properly read.
  - Altimeter Gridded netcdf files (generated by gridSatGlobal_Altimeter.py)
    that must be written into a list named satlist.txt. Each line contains
-   that must be written into a list named satlist.txt. Each line containing
    one satellite mission.
  - WAVEWATCHIII netcdf field outputs, also saved in a list ww3list.txt
  It is important to mention that the grid points of WAVEWATCHIII results
@@ -221,7 +220,6 @@
 	ncfile.history="Matchups of WAVEWATCHIII and AODN Altimeter data. Total of "+repr(size(ind))+" observations or pairs model/observation."
 	# create  dimensions. 2 Dimensions
 	ncfile.createDimension('index',ftime.shape[0])
-	ncfile.createDimension('time',ftime.shape[0])
 	ncfile.createDimension('satellite', sdname.shape[0] )
 	ncfile.createDimension('GlobalOceansSeas', ocnames.shape[0] )
 	ncfile.createDimension('HighSeasMarineZones', hsmznames.shape[0] )

validation/merrevalstats.py → validation/mvalstats.py

@@ -2,37 +2,36 @@
 # -*- coding: utf-8 -*-
 
 """
-merrevalstats.py
+mvalstats.py
 
 VERSION AND LAST UPDATE:
  v1.0  04/04/2022
+ v1.1  05/22/2022
 
 PURPOSE:
  Group of python functions for statistical analyses and validation. 
  Users can import as a standard python function, and use it accordingly:
  For example:
-  import merrevalstats
-  merrevalstats.metrics(model,obs,0.01,20,15)
+  import mvalstats
+  mvalstats.metrics(model,obs,0.01,20,15)
  Users can help() each function to obtain information about inputs/outputs
-  help(merrevalstats.metrics)
+  help(mvalstats.metrics)
 
 USAGE:
  functions
    smrstat
    metrics
- Explanation for each function is contained in the headers
+ Explanation of each function is contained in the headers
 
 OUTPUT:
  summary statistics values; and error metrics 
 
 DEPENDENCIES:
  See dependencies.py and the imports below.
- lmoments3 is a library/module not included in anaconda by default.
-  It can be installed with
-  pip install git+https://github.com/OpenHydrology/lmoments3.git
 
 AUTHOR and DATE:
  04/04/2022: Ricardo M. Campos, first version.
+ 05/22/2022: Ricardo M. Campos, removed L-moments from smrstat.
 
 PERSON OF CONTACT:
  Ricardo M Campos: ricardo.campos@noaa.gov
@@ -42,9 +41,8 @@
 from pylab import *
 from matplotlib.mlab import *
 import numpy as np
-import copy
 import scipy.stats
-import lmoments3 as lm
+# import lmoments3 as lm
 # pip install git+https://github.com/OpenHydrology/lmoments3.git
 import sys
 import warnings; warnings.filterwarnings("ignore")
@@ -56,16 +54,18 @@ def smrstat(*args):
 	Input: one array of interest
 	Output: mean, variance, skewness, kurtosis, min, max, percentile80, percentile90, percentile95, percentile99, percentile99.9
 	Example:
-	  import merrevalstats
-	  sresult = merrevalstats.smrstat(hs,0,20)
+	  import mvalstats
+	  sresult = mvalstats.smrstat(hs)
+	  sresult = mvalstats.smrstat(hs,0,20)
 	'''
+
 	vmin=-np.inf; vmax=np.inf
 	if len(args) == 1:
-		x=copy.copy(args[0])
+		x=np.array(args[0])
 	elif len(args) == 2:
-		x=copy.copy(args[0]); vmin=copy.copy(args[1])
+		x=np.array(args[0]); vmin=np.float(args[1])
 	elif len(args) == 3:
-		x=copy.copy(args[0]); vmin=copy.copy(args[1]); vmax=copy.copy(args[2])
+		x=np.array(args[0]); vmin=np.float(args[1]); vmax=np.float(args[2])
 	elif len(args) > 3:
 		sys.exit(' Too many inputs')
 
@@ -75,7 +75,8 @@ def smrstat(*args):
 	else:
 		sys.exit(' Array without valid numbers.')
 
-	ferr=np.zeros((14),'f')*np.nan
+	# ferr=np.zeros((14),'f')*np.nan
+	ferr=np.zeros((11),'f')*np.nan
 	ferr[0] = np.mean(x)
 	ferr[1] = np.var(x)
 	ferr[2] = scipy.stats.skew(x)
@@ -89,10 +90,10 @@ def smrstat(*args):
 	ferr[10] = np.percentile(x,99.9)
 	# Hosking & Wallis L-moment ratios
 	# pip install git+https://github.com/OpenHydrology/lmoments3.git
-	hwlm = lm.lmom_ratios(x, nmom=5)
-	ferr[11] = hwlm[1]/hwlm[0]
-	ferr[12] = hwlm[2]
-	ferr[13] = hwlm[3]
+	# hwlm = lm.lmom_ratios(x, nmom=5)
+	# ferr[11] = hwlm[1]/hwlm[0]
+	# ferr[12] = hwlm[2]
+	# ferr[13] = hwlm[3]
 
 	return ferr
 
@@ -107,22 +108,26 @@ def metrics(*args):
 		They must have the same size
 	Output: numpy array with shape equal to 8:
 		bias, RMSE, NBias, NRMSE, SCrmse, SI, HH, CC
+	Example:
+	  import mvalstats
+	  mresult = mvalstats.metrics(model_hs,obs_hs)
+	  mresult = mvalstats.metrics(model_hs,obs_hs,0,20)
 	'''
 
 	vmin=-np.inf; vmax=np.inf; maxdiff=np.inf
 	if len(args) < 2:
 		sys.exit(' Need two arrays with model and observations.')
 	elif len(args) == 2:
-		model=copy.copy(args[0]); obs=copy.copy(args[1])
+		model=np.array(args[0]); obs=np.array(args[1])
 	elif len(args) == 3:
-		model=copy.copy(args[0]); obs=copy.copy(args[1])
-		vmin=copy.copy(args[2])
+		model=np.array(args[0]); obs=np.array(args[1])
+		vmin=np.float(args[2])
 	elif len(args) == 4:
-		model=copy.copy(args[0]); obs=copy.copy(args[1])
-		vmin=copy.copy(args[2]); vmax=copy.copy(args[3])
+		model=np.array(args[0]); obs=np.array(args[1])
+		vmin=np.float(args[2]); vmax=np.float(args[3])
 	elif len(args) == 5:
-		model=copy.copy(args[0]); obs=copy.copy(args[1])
-		vmin=copy.copy(args[2]); vmax=copy.copy(args[3]); maxdiff=copy.copy(args[4]);
+		model=np.array(args[0]); obs=np.array(args[1])
+		vmin=np.float(args[2]); vmax=np.float(args[3]); maxdiff=np.float(args[4])
 	elif len(args) > 5:
 		sys.exit(' Too many inputs')
 

validation/organizeDistanceToCoast.py

@@ -38,6 +38,7 @@
 import matplotlib.pyplot as plt
 from mpl_toolkits.basemap import cm
 from netCDF4 import Dataset
+import warnings; warnings.filterwarnings("ignore")
 colormap = cm.GMT_polar
 palette = plt.cm.jet
 palette.set_bad('aqua', 10.0)

validation/prepGridMask.py

@@ -90,8 +90,7 @@
 import cartopy
 from matplotlib import ticker
 # import pickle
-import warnings
-warnings.filterwarnings("ignore")
+import warnings; warnings.filterwarnings("ignore")
 # netcdf format
 fnetcdf="NETCDF4"