Commit d3e21bbf authored by GROUAZEL's avatar GROUAZEL

modif for branch plot_spec_v3

parent b9dcfd9c
#from plot_spec_V2 import plot_spec_generic
from plot_spec_V2 import spec_plot
#from plot_spec_V2 import spec_plot
from plot_spec_V3 import spec_plot
import spectra_polar_plot_xarray #graphics
from plot_spec_V2 import getColorMap
import logging
import matplotlib.pyplot as plt
......@@ -173,8 +175,7 @@ class Spectrum:
"""
if not style:
style = init_spectrum_style(self.meta.type)
if self.meta.type[0] == "s":
#Test for imaginary part
if hasattr(self.spec_data,'k'):
......@@ -188,6 +189,15 @@ class Spectrum:
sp_log = None)
return None
def generate_v2(self,style=None,cut_off=50):
"""
version Dec 2019 to use polar_plot_xarray instead of the classic spec_plot()
:param style:
:param cut_off:
:return:
"""
def show(self,style=None,cut_off=50):
"""
Function that will call the generation of a spectrum figure then show it
......@@ -199,8 +209,7 @@ class Spectrum:
plt.close(myFig)
def write(self,outputfile,showresult=False,style=None,
overwrite=False,dpi=100):
overwrite=False,dpi=100):
if not(overwrite) and (os.path.isfile(outputfile)):
pass
logging.debug('file:'+outputfile+' already exist!')
......
......@@ -23,8 +23,8 @@ class AbstractSpectrumReader():
#-----------------------------------------------
class SpectrumReaderSat(AbstractSpectrumReader):
@staticmethod
def read_spectrum(ncf=None,mode='satfull',dataName='UNK'):
# @staticmethod
def read_spectrum(self,ncf=None,mode='satfull',dataName='UNK'):
"""
read_spectrum_sat: Function that return a spectrum object from a satelite file
ncf: Net CDF File already open
......@@ -48,16 +48,16 @@ class SpectrumReaderSat(AbstractSpectrumReader):
"""
#index is useless for now
myspec_data = SpectrumReaderSat.read_spec_data(ncf,mode)
mywsys = SpectrumReaderSat.read_wsys(ncf)
mywind = SpectrumReaderSat.read_wind(ncf)
myspec_data = self.read_spec_data(ncf,mode)
mywsys = self.read_wsys(ncf)
mywind = self.read_wind(ncf)
logging.debug('read_spectrum | mywind = %s',mywind)
mymeta = SpectrumReaderSat.read_meta(ncf,mode)
mymeta = self.read_meta(ncf,mode)
return spc.Spectrum(myspec_data,mywsys,mywind,mymeta,dataName=dataName)
@staticmethod
def read_spectrum_list(ncf=None):
# @staticmethod
def read_spectrum_list(self,ncf=None):
"""
read_spectrum_list_sat: Function that return an array with the spectrum objects from a satelite file
it will return the following:
......@@ -86,7 +86,8 @@ class SpectrumReaderSat(AbstractSpectrumReader):
fileToRead = ncf
for mode in modes:
myspectrum = SpectrumReaderSat.read_spectrum(ncf=fileToRead,mode=mode,dataName=dataName)
myspectrum = self.read_spectrum(ncf=fileToRead,mode=mode,dataName=dataName)
# myspectrum = SpectrumReaderSat.read_spectrum(ncf=fileToRead,mode=mode,dataName=dataName)
spectrumdict[mode] = myspectrum
if isinstance(ncf, basestring):
fileToRead.close()
......@@ -103,6 +104,11 @@ class SpectrumReaderSat(AbstractSpectrumReader):
wsys['dirad'] = ncf.variables['oswDirmet'][:][0][0]*np.pi/180.
wsys['lon'] = ncf.variables['oswLon'][:][0][0]
wsys['lat'] = ncf.variables['oswLat'][:][0][0]
if 'oswQualityFlagPartition' in ncf.variables:
wsys['qcpartition'] = ncf.variables['oswQualityFlagPartition'][:][0][0]
else:
wsys['qcpartition'] = np.ones((5,))*np.nan
# logging.debug('qc partition = %s',wsys['qcpartition'])
return wsys
@staticmethod
......@@ -120,13 +126,13 @@ class SpectrumReaderSat(AbstractSpectrumReader):
logging.debug('ESA direction given in convention FROM %s',dirwind)
dirwind = np.mod(dirwind+180,360)
logging.debug('direction given in convention TO %s',dirwind)
if False:
wind['U'] = ncf.variables['oswWindSpeed'][:][0][0]*np.cos(np.radians(dirwind)) #correction A. Grouazel Sept 2019
wind['V'] = ncf.variables['oswWindSpeed'][:][0][0]*np.sin(np.radians(dirwind))
else:
wind['U'] = ncf.variables['oswWindSpeed'][:][0][0]*np.cos(np.pi/2-np.radians(dirwind)) #correction A. Grouazel Sept 2019
wind['V'] = ncf.variables['oswWindSpeed'][:][0][0]*np.sin(np.pi/2-np.radians(dirwind))
wind['dir'] = dirwind
# if False:
# wind['U'] = ncf.variables['oswWindSpeed'][:][0][0]*np.cos(np.radians(dirwind)) #correction A. Grouazel Sept 2019
# wind['V'] = ncf.variables['oswWindSpeed'][:][0][0]*np.sin(np.radians(dirwind))
# else:
wind['U'] = ncf.variables['oswWindSpeed'][:][0][0]*np.cos(np.pi/2-np.radians(dirwind)) #correction A. Grouazel Sept 2019
wind['V'] = ncf.variables['oswWindSpeed'][:][0][0]*np.sin(np.pi/2-np.radians(dirwind))
wind['dir'] = dirwind
logging.debug('u = %s v = %s',wind['U'],wind['V'])
return wind
......@@ -165,14 +171,16 @@ class SpectrumReaderSat(AbstractSpectrumReader):
my_data['phi'] = ncf.variables['oswPhi'][:]
# print('mode1',mode,my_data['phi'],type(my_data['phi']))
if isinstance(my_data['phi'],np.ma.core.MaskedArray) or True:
if isinstance(my_data['phi'],np.ma.core.MaskedArray) or True: #je ne comprend pas a quoi sert ce test...
my_data['phi'] = np.ma.array(my_data['phi'])
my_data['phi'].mask = np.zeros(my_data['phi'].shape)
#special patch for oswK [agrouaze 11/06/2016]
# (values set to the fill value do not have the attribute fill=True thus they are not retrieved as masked values)
my_data['k'] = np.ma.masked_where(my_data['k']>=9.96920e+36, my_data['k'], copy=False)
if mode == "satfull":
my_data['sp'] = ncf.variables['oswPolSpec'][:][0][0]
tmptmp =ncf.variables['oswPolSpec'][:][0][0]
my_data['sp'] = np.ma.masked_where(tmptmp >= 10000,tmptmp,copy = False)
print('spectre ',my_data['sp'].shape)
if not isinstance(my_data['k'].mask,np.ndarray):
my_data['k'].mask = np.zeros(my_data['k'].shape)
......
# coding: utf-8
import abc
import numpy as np
import netCDF4 as nc
......@@ -8,13 +9,16 @@ import Spectrum as spc
import logging
import spectrumUtil as spUtil
import os
from ocean_wave_spectrum_partitioning_IPF_proto import partition_spectrum #quality
from reconstruct_2D_partition_label_matrix import reconstruct_ww3_partition_mask_from_partitioning #quality
from convert_partitioning_output_2_wsys_polar_spectral_plot import convert_partitioning_output_2_wsys
from get_ww3grid_fromS1measu_L2_fullpath import get_path_ww3grid,get_closest_gridpoint #colocation
class AbstractSpectrumReader():
__metaclass__ = abc.ABCMeta
@abc.abstractmethod
def read_spectrum(ncfile,mode="satfull",index=None):
def read_spectrum(self,ncfile,mode="satfull",index=None):
pass
#-----------------------------------------------
......@@ -22,8 +26,9 @@ class AbstractSpectrumReader():
#-----------------------------------------------
class SpectrumReaderWw3(AbstractSpectrumReader):
@staticmethod
def read_spectrum(ncf=None,mode='ww3',dataName='UNK'):
# @staticmethod
def read_spectrum(self,ncf=None,mode='ww3',dataName='UNK',path=None):
"""
read_spectrum_ww3: Function that return a spectrum object from a WW3 prevision file
ncf: Net CDF File already open
......@@ -42,22 +47,43 @@ class SpectrumReaderWw3(AbstractSpectrumReader):
myspectrum = read_spectrum_ww3(path='path/to/ncfile.nc')
"""
myspec_data = SpectrumReaderWw3.read_spec_data(ncf)
mywsys = SpectrumReaderWw3.read_wsys(ncf)
mywind = SpectrumReaderWw3.read_wind(ncf)
mymeta = SpectrumReaderWw3.read_meta(ncf,mode)
if path is not None:
self.fullpath = path
else:
self.fullpath = None
myspec_data = self.read_spec_data(ncf)
mywsys = self.read_wsys(ncf)
self.mywsys = mywsys
mywind = self.read_wind(ncf)
mymeta = self.read_meta(ncf,mode)
# myspec_data = SpectrumReaderWw3.read_spec_data(ncf)
# mywsys = SpectrumReaderWw3.read_wsys(ncf)
# mywind = SpectrumReaderWw3.read_wind(ncf)
# mymeta = SpectrumReaderWw3.read_meta(ncf,mode)
return spc.Spectrum(myspec_data,mywsys,mywind,mymeta,dataName=dataName)
@staticmethod
def read_wsys(ncf):
# @staticmethod
def read_wsys(self,ncf):
wsys = {}
wsys['lon'] = ncf.variables['lon'][:]
wsys['lon'] = float(wsys['lon'])
wsys['lat'] = ncf.variables['lat'][:]
wsys['lat'] = float(wsys['lat'])
# wsys['lon'] = ncf.variables['lon'][:]
# logging.debug('SpectrumReaderWw3 | wsys["lon"] = %s',wsys['lon'])
wsys['lon'] = self.lon
# wsys['lat'] = ncf.variables['lat'][:]
wsys['lat'] = self.lat
wsys['hs'] = self.tmp_wsys['hs']
wsys['wl'] = self.tmp_wsys['wl']
wsys['dirad'] = self.tmp_wsys['dirad']
wsys['hsgrid'] = self.hsgrid
return wsys
@staticmethod
def read_wind(ncf):
# def do_partitioning(ncf):
# return complet_partitions_sorted,wave_integrated_values
# @staticmethod
def read_wind(self,ncf):
"""
since Mickael Accensi wind in colocated WW3 spectra a in FROM convention 'example: +90° mean coming from east'
#which is the meteorological convention also used in the new IPF2.90 since 13march2018
......@@ -73,28 +99,55 @@ class SpectrumReaderWw3(AbstractSpectrumReader):
# wind['V'] = -wind['V']
if True:#patch agrouaze Sept 2019: (90-val+180)%360 correction sans explication qui vient de tests sur http://grougrou1.ifremer.fr:8888/notebooks/workspace_leste/cfosat_perso/verif_windir_computation.ipynb#
# wdir = (90-np.arctan2(wind['U'],wind['V']) * 180 / np.pi + 180.) % 360
wdir = (np.arctan2(wind['V'],wind['U']) * 180 / np.pi + 180.) % 360 #trop bien jai trouvé que en inversant V et U comme dans la doc numpy on retombe bien sur la bonne direction
wdir = (np.arctan2(wind['V'],wind['U']) * 180 / np.pi + 180.) % 360 #en inversant V et U comme dans la doc numpy on retombe bien sur la bonne direction
#Il faut mettre le V d abord (ie la composent meridional) et en second le U (ie composante zonale)
wind['V'] = speed*np.cos(np.radians(wdir))
wind['U'] = speed*np.sin(np.radians(wdir))
# print('plop')
else:
wdir = (np.arctan2(wind['U'],wind['V']) * 180 / np.pi) % 360
wind['dir'] = wdir
return wind
@staticmethod
def read_meta(ncf,mode):
# @staticmethod
def read_meta(self,ncf,mode):
meta={}
meta['type'] = mode
return meta
@staticmethod
def read_spec_data(ncf):
# @staticmethod
def read_spec_data(self,ncf):
my_data = {}
my_data['k'] = ncf.variables['k'][:]
my_data['phi'] = ncf.variables['phi'][:]
my_data['sp'] = ncf.variables['polSpec'][:].squeeze()
self.lon = float(ncf.variables['lon'][:])
self.lat = float(ncf.variables['lat'][:])
area = ncf.variables['area'][:]
my_data['sp'] = ncf.variables['polSpec'][:].squeeze() #*area added by agrouaze Nov 2019 to have reasonable Hs at the partitioning step
#partitioning
#valeur utilisees pour les cross assignment des bulletins mensuels Dec 2019
noise_thr = 2.
noise_thr = 0.1
merge_thr = 4.
discard_thr = 0.
#test avec une autre valeur de discard
merge_thr = 2
discard_thr = 8
# discard_thr = 0.022
wavesys = partition_spectrum(my_data['sp'],my_data['k'], my_data['phi'], npartitions=5, twod=True, symphi=False,
klow=None, khigh=None,
noise_thr=noise_thr, merge_thr=merge_thr,
discard_thr=discard_thr, polar=True, alfah=None, merge_close=True,
rel_peak_high=False, spec_std0=None, rem=None,area=area)
complet_partitions_sorted = reconstruct_ww3_partition_mask_from_partitioning(my_data['sp']*area,wavesys)
# complet_partitions_sorted = partition_wv_xassignment.reconstruct_ww3_partition_mask_from_partitioning(complet_matrice_sorted,ww3_part)
wave_integrated_values = convert_partitioning_output_2_wsys(wavesys)
#end partitioning
my_data['part'] = complet_partitions_sorted
t = ncf.variables['time'][:]
u = ncf.variables['time'].units
my_data['start_time'] = nc.num2date(t,u)
......@@ -104,7 +157,21 @@ class SpectrumReaderWw3(AbstractSpectrumReader):
u = ncf.variables['sar_time'].units
my_data['sar_time'] = nc.num2date(t,u)
my_data['sar_time'] = my_data['sar_time'][0]
# my_data['tmp_wsys'] = wave_integrated_values
#add hs from closest ww3 grid point
if self.fullpath is not None:
res_file = get_path_ww3grid(measuL2WV_fullpath = self.fullpath)
handler_ww3gridfile = Dataset(res_file)
idlat,idlon = get_closest_gridpoint(self.lon,self.lat,handler_ww3gridfile)
my_data['hsgrid'] = handler_ww3gridfile.variables['hs'][0,idlat,idlon]
handler_ww3gridfile.close()
else:
my_data['hsgrid'] = np.nan
logging.debug('ww3 hs from grid = %sm',my_data['hsgrid'])
self.tmp_wsys = wave_integrated_values
self.hsgrid = my_data['hsgrid']
return my_data
#-----------------------------------------------
# End of Reader for WW3 File
......
......@@ -12,22 +12,24 @@ from SpectrumReaderBuoys import SpectrumReaderBuoyGlobwave
from SpectrumReaderSentinel import SpectrumReaderSat
from SpectrumReaderWW3 import SpectrumReaderWw3
class AbstractSpectrumReader():
__metaclass__ = abc.ABCMeta
@abc.abstractmethod
def read_spectrum(ncfile,mode="satfull",index=None):
# @abc.abstractmethod
def read_spectrum(self,ncfile,mode="satfull",index=None):
pass
instww3 = SpectrumReaderWw3()
insts1 = SpectrumReaderSat()
instbuoy = SpectrumReaderBuoyGlobwave()
readerdictionary ={
"satreal":SpectrumReaderSat.read_spectrum,
"satfull":SpectrumReaderSat.read_spectrum,
"satimaginary":SpectrumReaderSat.read_spectrum,
"ww3":SpectrumReaderWw3.read_spectrum,
"buoyglobwavendbc":SpectrumReaderBuoyGlobwave.read_spectrum,
"buoyglobwavecdip":SpectrumReaderBuoyGlobwave.read_spectrum
"satreal":insts1.read_spectrum,
"satfull":insts1.read_spectrum,
"satimaginary":insts1.read_spectrum,
"ww3":instww3.read_spectrum,
"buoyglobwavendbc":instbuoy.read_spectrum,
"buoyglobwavecdip":instbuoy.read_spectrum
}
......@@ -36,25 +38,30 @@ readerdictionary ={
# Generic Reader that call the good reader
#-----------------------------------------------
class SpectrumReader(AbstractSpectrumReader):
@staticmethod
def read_spectrum(ncfile,mode=None,index=None):
def __init__(self):
self.fullpath = None
# @staticmethod
def read_spectrum(self,ncfile,mode=None,index=None):
dataName='UNK'
if isinstance(ncfile, basestring):
dataName = SpectrumReader.selectDataName(ncfile)
if isinstance(ncfile, str):
# dataName = SpectrumReader.selectDataName(ncfile)
self.fullpath = ncfile
dataName = self.selectDataName(ncfile)
fileToRead = Dataset(ncfile, 'r')
if not mode:
mode = SpectrumReader.selectmodebyname(ncfile)
# mode = SpectrumReader.selectmodebyname(ncfile)
mode = self.selectmodebyname(ncfile)
else:
fileToRead = ncfile
SpectrumReader.selectmodebydata(ncfile)
# SpectrumReader.selectmodebydata(ncfile)
self.selectmodebydata(ncfile)
if index:
myspectrum = readerdictionary[mode](fileToRead,index=index,mode=mode,dataName=dataName)
else:
myspectrum = readerdictionary[mode](fileToRead,mode=mode,dataName=dataName)
if isinstance(ncfile, basestring):
if isinstance(ncfile, str):
fileToRead.close()
#myspectrum.dataName=dataName
......@@ -62,8 +69,8 @@ class SpectrumReader(AbstractSpectrumReader):
return myspectrum
@staticmethod
def selectDataName(name):
# @staticmethod
def selectDataName(self,name):
if 'globwave' in name and 'cdip' in name:
return 'buoyglobwavecdip'
elif 'globwave' in name and 'ndbc' in name:
......@@ -77,8 +84,8 @@ class SpectrumReader(AbstractSpectrumReader):
@staticmethod
def selectmodebyname(name):
# @staticmethod
def selectmodebyname(self,name):
if 'globwave' in name and 'cdip' in name:
return 'buoyglobwavecdip'
elif 'globwave' in name and 'ndbc' in name:
......@@ -88,8 +95,8 @@ class SpectrumReader(AbstractSpectrumReader):
elif 's1a' in name or 's1b' in name:
print('ambiguous situation will return the full spectrum mode:satfull')
return 'satfull'
@staticmethod
def selectmodebydata(data):
# @staticmethod
def selectmodebydata(self,data):
print(data.title)
#-----------------------------------------------
# End of Generic Reader
......
Markdown is supported
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