Commit ac5ca3c4 authored by GROUAZEL's avatar GROUAZEL

fix issue of wind direction for S1 and WW3

parent 5dd1cb03
#from plot_spec_V2 import plot_spec_generic
from plot_spec_V2 import spec_plot
from plot_spec_V2 import getColorMap
import logging
import matplotlib.pyplot as plt
import matplotlib as mpl
import os
......@@ -167,7 +167,7 @@ class Spectrum:
else:
self.title = title
def generate(self,style=None):
def generate(self,style=None,cut_off=50):
"""
Function that will call the generation of a spectrum figure then return it
"""
......@@ -179,20 +179,20 @@ class Spectrum:
#Test for imaginary part
if hasattr(self.spec_data,'k'):
if isinstance(self.spec_data.k,np.ma.core.MaskedArray):
return spec_plot(self,style,cut_off = 100,
return spec_plot(self,style,cut_off = cut_off,
sp_log = None)
else:
plt.annotate('Spectrum cannot be displayed because of unusable values',xy=(0.375, 0.5),color=(1., .4, 0.))
else:
return spec_plot(self,style, cut_off = 100,
return spec_plot(self,style, cut_off = cut_off,
sp_log = None)
return None
def show(self,style=None):
def show(self,style=None,cut_off=50):
"""
Function that will call the generation of a spectrum figure then show it
"""
myFig = self.generate(style=style)
myFig = self.generate(style=style,cut_off=cut_off)
#myFig = plot_spec_generic(self,mycolormap=mycolormap,dateFormat=dateFormat)
plt.show()
if myFig:
......@@ -203,7 +203,7 @@ class Spectrum:
if not(overwrite) and (os.path.isfile(outputfile)):
pass
print('file:'+outputfile+' already exist!')
logging.debug('file:'+outputfile+' already exist!')
else:
myFig = self.generate(style=style)
if 'eps' in outputfile:
......
......@@ -51,6 +51,7 @@ class SpectrumReaderSat(AbstractSpectrumReader):
myspec_data = SpectrumReaderSat.read_spec_data(ncf,mode)
mywsys = SpectrumReaderSat.read_wsys(ncf)
mywind = SpectrumReaderSat.read_wind(ncf)
logging.debug('read_spectrum | mywind = %s',mywind)
mymeta = SpectrumReaderSat.read_meta(ncf,mode)
return spc.Spectrum(myspec_data,mywsys,mywind,mymeta,dataName=dataName)
......@@ -108,10 +109,25 @@ class SpectrumReaderSat(AbstractSpectrumReader):
def read_wind(ncf):
"""
read_wind_sat: Function that return a wind dictionary from a sat file
oswWindDirection:long_name = "SAR Wind direction (meteorological convention=clockwise direction from where the wind comes respect to North)"
note: avant IPF 2.91 la convention etait erronee
"""
wind={}
wind['U'] = ncf.variables['oswWindSpeed'][:][0][0]*np.cos(90- ncf.variables['oswWindDirection'][:][0][0]*np.pi/180.)
wind['V'] = ncf.variables['oswWindSpeed'][:][0][0]*np.sin(90- ncf.variables['oswWindDirection'][:][0][0]*np.pi/180.)
# wind['U'] = ncf.variables['oswWindSpeed'][:][0][0]*np.cos(90- ncf.variables['oswWindDirection'][:][0][0]*np.pi/180.) #original version from MLB ?
# wind['V'] = ncf.variables['oswWindSpeed'][:][0][0]*np.sin(90- ncf.variables['oswWindDirection'][:][0][0]*np.pi/180.)
# oswWindDirection:long_name = "SAR Wind direction (meteorological convention=clockwise direction from where the wind comes respect to North)"
dirwind = ncf.variables['oswWindDirection'][:][0][0]
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
logging.debug('u = %s v = %s',wind['U'],wind['V'])
return wind
......@@ -148,7 +164,8 @@ class SpectrumReaderSat(AbstractSpectrumReader):
my_data['k'] = ncf.variables['oswK'][:]
my_data['phi'] = ncf.variables['oswPhi'][:]
if not isinstance(my_data['phi'],np.ma.core.MaskedArray):
print('mode1',mode,my_data['phi'],type(my_data['phi']))
if isinstance(my_data['phi'],np.ma.core.MaskedArray) or True:
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]
......@@ -156,7 +173,6 @@ class SpectrumReaderSat(AbstractSpectrumReader):
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]
if not isinstance(my_data['k'].mask,np.ndarray):
my_data['k'].mask = np.zeros(my_data['k'].shape)
......@@ -185,6 +201,8 @@ class SpectrumReaderSat(AbstractSpectrumReader):
my_data['sp'] = my_data['sp'][:,kmask]
#recuperation sar_start_time
my_data['start_time'] = parser.parse(ncf.getncattr('firstMeasurementTime'))
else:
print('rien')
return my_data
#-----------------------------------------------
# End of Reader for Sat File
......
# coding: utf-8
import abc
import numpy as np
import netCDF4 as nc
......@@ -66,6 +67,8 @@ class SpectrumReaderWw3(AbstractSpectrumReader):
wind['U'] = wind['U'][0]
wind['V'] = ncf.variables['wind_speed_model_v'][:]
wind['V'] = wind['V'][0]
wdir = (np.arctan2(wind['U'],wind['V']) * 180 / np.pi) % 360
wind['dir'] = wdir
return wind
@staticmethod
def read_meta(ncf,mode):
......
......@@ -73,7 +73,7 @@ def init_sp(k,phi,spec):
#sp = np.ma.array(sp)
else:
logging.debug('on fait pas de manip sur le spectre')
print 'on fait pas de manip sur le spectre'
print('on fait pas de manip sur le spectre')
sp = spec
# converting the direction to rad
......@@ -93,7 +93,10 @@ def init_ax(ax):
ax.radii_angle = 90. # position of label name on radial axis
ax.theta_angles = np.arange(0, 360, 90)
ax.theta_labels = ['N', 'E', 'S', 'W']
ax.set_thetagrids(angles=ax.theta_angles, labels=ax.theta_labels,frac=1.02)
try:
ax.set_thetagrids(angles=ax.theta_angles, labels=ax.theta_labels,frac=1.02)
except:
ax.set_thetagrids(angles=ax.theta_angles, labels=ax.theta_labels)
ax.set_axis_off
ax.set_frame_on(False)
ax.xaxis.set_visible(True)
......@@ -382,22 +385,56 @@ def add_basic_elements_to_plot(ax_plot,title,fs_title,circleradius,freq):
#Function that add a wind Vector to ax_plot
def add_wind_to_plot(myspectrum,ax_plot,cut_off):
print('myspectrum',myspectrum)
print('add_wind_to_plot')
scriptpath = os.path.realpath(__file__)
print("Script path is : " + scriptpath)
U = None
V = None
#Test to ignore empty array
if myspectrum.wind and hasattr(myspectrum.wind,"U")and hasattr(myspectrum.wind,"V"):
U = myspectrum.wind.U #zonal
V = myspectrum.wind.V #meridional
logging.debug('get U and V from myspectrum')
else:
logging.debug('myspectrum.wind = %s',myspectrum.wind)
# plot the wind vector
if U and V:
wspd = np.sqrt(U**2+V**2)
ax_plot.quiver(np.arctan2(U,V),0.5*2*np.pi/cut_off,
U/wspd,V/wspd,
scale=5,width=0.015,color='red')
ax_plot.annotate('U10 = {:3.2f}'.format(wspd)+ ' m/s',
if False: #ca marche pour mon cas mais je suis quasiment sur que ca va faire nimp dans les autres cas
wspd = np.sqrt(U**2+V**2)
x = np.arctan2(-U,-V) #position in radians
y = 0.5*2*np.pi/cut_off #distance from the center of the plot
logging.debug('x=%s y=%s',x,y)
logging.debug('x in deg = %s',np.degrees(x))
ax_plot.quiver(x,y,
-U/wspd,-V/wspd,
scale=5,width=0.015,color='red')
else:
#a condition que les U et V aient ete calculer comme ca:U = 10*np.cos(np.pi/2-np.radians(dirwind))
wspd = np.sqrt(U**2+V**2)
x = np.arctan2(U,V) #position in radians
y = 0.5*2*np.pi/cut_off #distance from the center of the plot
logging.debug('x=%s y=%s',x,y)
logging.debug('x in deg = %s',np.degrees(x))
ax_plot.quiver(x,y,
U/wspd,V/wspd,
scale=5,width=0.015,color='red')
# ax_plot.annotate(r'$U10_TO$ = {:3.2f}'.format(wspd)+ ' m/s',
if 'dir' in dir(myspectrum.wind):
dire = myspectrum.wind.dir
else:
dire = np.nan
try:
ax_plot.annotate(r'$U10_{TO}$ = %3.2f'%(wspd)+ ' m/s\ndir=%3.1f$^o$'%(myspectrum.wind.dir),
xy=(np.arctan2(U,V),0.9*2*np.pi/cut_off),color='red',zorder=1000000)
except:
ax_plot.annotate(r'$U10_{TO}$ = %3.2f'%(wspd)+ ' m/s',
xy=(np.arctan2(U,V),0.9*2*np.pi/cut_off),color='red',zorder=1000000)
logging.debug('quiver added')
else:
logging.debug('quiver not added')
return ax_plot
#Draw the spectrum and the contours
......@@ -447,7 +484,7 @@ def add_spectrum_to_plot(fig,myspectrum,vmin,vmax,colormap,cb_label,nrow,ncolumn
#adding the colorbar
fig = add_colorbar(fig,cb_label,vmax = vmax,vmin = vmin,colormap= colormap,nrow=nrow,ncolumn=ncolumn)
else:
print 'warning on values spectra'
print('warning on values spectra')
return fig,part_u
# Function that will generate a figure for a spectrum with all elements.
......
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