Commit 748c1613 authored by GROUAZEL's avatar GROUAZEL
Browse files

libs for cross spectra computation from SLC

parent 31f6d181
......@@ -21,7 +21,7 @@ import xsarsea
project = 'xsarsea'
copyright = '2020, Ifremer LOPS/SIAM'
author = 'Olivier Archer, Alexandre Levieux'
version = xsarsea.xsarsea.__version__
version = xsarsea.detrend.__version__
# -- General configuration ---------------------------------------------------
......
from .xsarsea import *
# from .streaks import wind_streaks
#!/usr/bin/env python
import numpy as np
import numexpr as ne
import xarray as xr
import time
import sys
......@@ -27,6 +28,10 @@ def get_imagette_indice(onetiff,wv_slc_meta):
logging.info('indice in meta sentinel1 gdal driver is : %s',good_indice)
return good_indice
from numba import jit
@jit(nopython=True,parallel=True)
def supermean(a):
return np.mean(a)
def read_annotations(annotation_file_xml):
"""
......@@ -49,33 +54,46 @@ def read_annotations(annotation_file_xml):
#meta_dict = {PIL.TiffTags.TAGS[key] : img.tag[key] for key in img.tag_v2.keys()}
return rangeSpacing,azimuthSpacing,mean_inc
def read_slc(onetiff,slice_subdomain=None):
def read_slc(onetiff,slice_subdomain=None,resolution=None,resampling=None):
"""
:param one_wv: SAFE (str) full path without the last /
:param slice_subdomain : slice objet to define the sub part of imagette to perform cross spectra
:param resolution: dict for instance for 200m {'atrack' : int(np.round(200 / sar_meta.pixel_atrack_m)), 'xtrack': int(np.round(200 / sar_meta.pixel_xtrack_m))}
:param rasterio.enums.Resampling.rms for instance
:return:
"""
logging.info('tiff: %s',onetiff)
logging.info('resampling : %s',resampling)
safepath = os.path.dirname(os.path.dirname(onetiff))
logging.info('safepath: %s',safepath)
tmp = xsar.get_test_file(safepath)
logging.info('tmp: %s',tmp)
wv_slc_meta = xsar.SentinelMeta(tmp)
wv_slc_meta = xsar.Sentinel1Meta(tmp)
logging.debug('wv_slc_meta %s',wv_slc_meta)
good_indice = get_imagette_indice(onetiff,wv_slc_meta)
slc = xsar.open_dataset(wv_slc_meta.subdatasets[good_indice])
if resolution is not None:
if resolution['atrack']==1 and resolution['xtrack']==1:
#June 2021, a patch because currently resolution 1:1 for image and rasterio returns an error
slc = xsar.open_dataset(wv_slc_meta.subdatasets[good_indice],resolution=None,resampling=None)
else:
slc = xsar.open_dataset(wv_slc_meta.subdatasets[good_indice],resolution=resolution,resampling=resampling)
else:
slc = xsar.open_dataset(wv_slc_meta.subdatasets[good_indice])
logging.info('max xtrack val : %s',slc['xtrack'].values[-1])
slc = slc.rename({'atrack' : 'azimuth','xtrack' : 'range'})
annotation_file_xml = onetiff.replace('measurement',"annotation").replace('.tiff','.xml')
rangeSpacing,azimuthSpacing,mean_inc = read_annotations(annotation_file_xml)
#rangeSpacing,azimuthSpacing,mean_inc = read_annotations(annotation_file_xml) # I prefere to use the result from xsar since this atribute is update with a specific resolution given
rangeSpacing = slc.attrs['pixel_xtrack_m']
azimuthSpacing= slc.attrs['pixel_atrack_m']
mean_inc = np.mean(slc['incidence']).values
platform_heading = slc.attrs['platform_heading']
logging.info('sland rangeSpacing %s',rangeSpacing)
logging.debug('slant rangeSpacing %s',rangeSpacing)
#changement des coordintes range and azimuth like Nouguier
raxis = np.arange(len(slc['range']))
aaxis = np.arange(len(slc['azimuth']))
logging.info('mean incidence : %s',mean_inc)
logging.info('factor to turn inout ground range: %s',rangeSpacing / np.sin(np.radians(mean_inc)))
logging.debug('mean incidence : %s',mean_inc)
logging.debug('factor to turn into ground range: %s',rangeSpacing / np.sin(np.radians(mean_inc)))
raxis = raxis * rangeSpacing / np.sin(np.radians(mean_inc))
aaxis = aaxis * azimuthSpacing
logging.info('range coords before : %s',slc['range'].values)
......@@ -154,7 +172,7 @@ def ground_regularization(ds, *,method='nearest', **kwargs):
def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0., look_window=None, range_spacing=None,
welsh_window='hanning', nperseg={'range':512, 'azimuth':512}, noverlap={'range':256, 'azimuth':256},
spacing_tol=1e-3,debug_plot=False,merge_along_tau=True, **kwargs):
spacing_tol=1e-3,debug_plot=False,merge_along_tau=True,return_periodoXspec=True, **kwargs):
"""
Compute SAR cross spectrum using a 2D Welch method. Looks are centered on the mean Doppler frequency
If ds contains only one cycle, spectrum wavenumbers are added as coordinates in returned DataSet, othrewise, they are passed as variables (k_range, k_azimuth).
......@@ -170,6 +188,7 @@ def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0.
look_window (xarray): window used in look processing
ds (xarray): R3S simulation or SLC
merge_along_tau (bool) : default is True -> final dimension should be like this (0tau: 48, 1tau: 32, 2tau: 16, freq_azimuth: 1024, freq_range: 1024, pol: 1) if False -> merge tau per sub region
return_periodoXspec (bool): default is True -> list of the X spectra (Im,Re 0,1,2 tau) returned for each sub images
kwargs (dict): other arguments passed to cowelch2d()
Returns:
......@@ -178,7 +197,7 @@ def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0.
if N_look<1:raise ValueError('N_look must be greater than 0')
if np.abs(look_width)>=1:raise ValueError('look_width must be in [0,1] range')
logging.debug('mean of modulation on the wole image %s',np.mean(slc['modulation'].values))
#logging.debug('mean of modulation on the wole image %s',np.mean(slc['modulation'].values))
with xr.set_options(keep_attrs=True):
gslc = ground_regularization(slc, range_spacing = range_spacing, **kwargs)
......@@ -219,8 +238,8 @@ def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0.
frange = np.fft.fftfreq(nperseg['range'], dir[0]/(2*np.pi))
fazimuth = np.fft.fftfreq(nperseg['azimuth'], dia[0]/(2*np.pi))
xspecs = list()
#variable_name = 'digital_number'
variable_name = 'modulation'
variable_name = 'digital_number'
#variable_name = 'modulation'
logging.debug('variable_name : %s',variable_name)
cpt=0
splitting_image ={}
......@@ -244,7 +263,7 @@ def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0.
logging.debug('axis_id_range = %s',axis_id_range)
logging.debug('sub au debut : %s',np.mean(sub[variable_name].values))
splitting_image[cpt] = sliceors
logging.debug('nb not fininte in s0 %s/%s ',(np.isfinite(sub[variable_name].values)==False).sum(),
logging.debug('nb not finite in s0 %s/%s ',(np.isfinite(sub[variable_name].values)==False).sum(),
sub[variable_name].values.size)
#s
if debug_plot:
......@@ -264,9 +283,19 @@ def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0.
#sub = sub[variable_name].values
logging.debug('sub after detrend: %s max: %s',np.mean(sub[variable_name].values),np.amax(sub[variable_name].values))
if cpt%20==0:
logging.info('sub %s %s %s/%s elapsed %1.1f sec',type(sub),sub[variable_name].shape,cpt,len(indicesx)*len(indicesy),time.time()-t0)
sub = sub[variable_name]/np.abs(sub[variable_name]).mean(dim=['range', 'azimuth']) # proposition Nouguier 17 may to do normalization before /detrend
#sub = sub[variable_name]/np.abs(sub[variable_name]).mean(dim=['range', 'azimuth']) # proposition Nouguier 17 may to do normalization before /detrend
tmpmean = supermean(np.abs(sub[variable_name]).values) # agrouaze test with numba to speed up
sub = sub[variable_name] /tmpmean
# logging.info('test numexpr : %s %s',type(sub[variable_name].values),sub[variable_name].values.dtype)
# tmpabs_re = ne.evaluate('abs(aa)',local_dict={'aa':sub[variable_name].values}).real
# tmpabs_1j = ne.evaluate('abs(aa)',local_dict={'aa':sub[variable_name].values}).imag
# tmpabs = tmpabs_re+1j*tmpabs_1j
# logging.info('tmpabs : %s %s %s',tmpabs.dtype,type(tmpabs),tmpabs.shape)
# tmpmean = ne.evaluate('sum(bb)',local_dict={'bb':tmpabs})/tmpabs.size
# sub = sub[variable_name]/tmpmean
# sub = sub*win # already commented in the original Nouguier code
if debug_plot and False:
plt.figure()
......@@ -313,9 +342,12 @@ def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0.
subxr = xr.DataArray(sub,dims=['pol','azimuth','range'])
else:
subxr = xr.DataArray(sub,dims=['azimuth','range'])
t1 = time.time()
tmplooks = compute_looks(subxr, N_look=N_look, look_width=look_width, look_overlap=look_overlap,
look_window=look_window, plot=plot,debug_plot=debug_plot)
if cpt % 20 == 0:
logging.info('sub %s %s %s/%s time compute_looks %1.3f sec time total to treat the sub image %1.3f sec',
type(sub),subxr.shape,cpt,len(indicesx)*len(indicesy),time.time()-t1,time.time()-t0)
#ajout agrouaze pour ajouter une coordonnee de portion dimage
for taux in tmplooks.keys():
for lookid,look in enumerate(tmplooks[taux]):
......@@ -347,7 +379,6 @@ def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0.
allspecs = list()
logging.info('N looks: %s',N_look)
logging.debug('xspecs : %s',xspecs)
#add allspecs per subdomain agrouaze
allspecs_per_sub_domain = {}
for subdomain in range(cpt):
......@@ -362,11 +393,12 @@ def compute_SAR_cross_spectrum(slc,*, N_look=3, look_width=0.25, look_overlap=0.
l.append(item.drop_vars(['freq_range','freq_azimuth']))
l = xr.concat(l,dim=str(tau) + 'tau').rename('cross-spectrum_' + str(tau) + 'tau')
all_spectra_one_domain.append(l)
all_spectra_one_domain = xr.merge(all_spectra_one_domain,compat='override')
all_spectra_one_domain = all_spectra_one_domain.assign_coords(freq_range=np.fft.fftshift(frange))
all_spectra_one_domain = all_spectra_one_domain.assign_coords(freq_azimuth=np.fft.fftshift(fazimuth))
all_spectra_one_domain = all_spectra_one_domain.rename(freq_range='kx',freq_azimuth='ky')
allspecs_per_sub_domain[subdomain] = all_spectra_one_domain
if return_periodoXspec:
all_spectra_one_domain = xr.merge(all_spectra_one_domain,compat='override')
all_spectra_one_domain = all_spectra_one_domain.assign_coords(freq_range=np.fft.fftshift(frange))
all_spectra_one_domain = all_spectra_one_domain.assign_coords(freq_azimuth=np.fft.fftshift(fazimuth))
all_spectra_one_domain = all_spectra_one_domain.rename(freq_range='kx',freq_azimuth='ky')
allspecs_per_sub_domain[subdomain] = all_spectra_one_domain
#add the limits of the sub images
limits_sub_images = xarray.Dataset()
......@@ -401,7 +433,7 @@ def compute_looks(gslc,*, N_look=3, look_width=0.25, look_overlap=0., look_windo
nlook = int(look_width*Np) # number of pulse in a look
noverlap = int(np.rint(look_overlap*look_width*Np)) # This noverlap is different from the noverlap of welch_kwargs
dim_range = 'range'
mydop = xrft.dft(gslc, dim=['azimuth'], detrend=None, window=None, shift=False,true_amplitude=true_ampl,true_phase=True)
mydop = xrft.dft(gslc, dim=['azimuth'], detrend=None, window=None, shift=False,true_amplitude=true_ampl,true_phase=False)
if debug_plot:
vmin = np.amin(mydop.squeeze().real.values)
vmax = np.amax(mydop.squeeze().real.values)
......@@ -443,9 +475,10 @@ def compute_looks(gslc,*, N_look=3, look_width=0.25, look_overlap=0., look_windo
mywin[{'freq_azimuth':slice(ind, ind+nlook)}]=win
if plot:mywin.plot()
# mywin = mywin.reindex_like(mydop) # This line is useless because it is automatically done when multiplying with mydop (xarray capabilities !)
look = xrft.idft(mydop*mywin, true_phase=True, dim=['freq_azimuth'], detrend=None, window=False, #true hase= True is test from agrouaze
look = xrft.idft(mydop*mywin, true_phase=False, dim=['freq_azimuth'], detrend=None, window=False, #true hase= True is test from agrouaze
shift=True,true_amplitude=true_ampl)
logging.debug('look mean values :%s',np.mean(look.values))
#logging.info('look. %s',look.freq_azimuth.spacing)
if debug_plot:
plt.figure()
plt.title('look #%s after idft (retour dans le domain spatial)'%ind)
......
......@@ -2,7 +2,7 @@ import logging
from pkg_resources import get_distribution
import numpy as np
__version__ = get_distribution("xsarsea").version
#__version__ = get_distribution("xsarsea").version
logger = logging.getLogger("xsarsea")
logger.addHandler(logging.NullHandler())
......
import sys
import time
import os
import logging
from matplotlib import pyplot as plt
from matplotlib import colors as mcolors
cmap = mcolors.LinearSegmentedColormap.from_list("", ["white","violet","mediumpurple","cyan","springgreen","yellow","red"])
import numpy as np
sys.path.append('/home1/datahome/agrouaze/git/xsarseafork/src/')
sys.path.append('/home1/datahome/agrouaze/git/xsarseafork/src/xsarsea/')
import xsarsea.conversion_polar_cartesian
import spectrum_clockwise_to_trigo
import display_xspectra_grid
import spectrum_rotation
reference_oswK_1145m_60pts = np.array([0.005235988, 0.00557381, 0.005933429, 0.00631625, 0.00672377,
0.007157583, 0.007619386, 0.008110984, 0.008634299, 0.009191379,
0.0097844, 0.01041568, 0.0110877, 0.01180307, 0.01256459, 0.01337525,
0.01423822, 0.01515686, 0.01613477, 0.01717577, 0.01828394, 0.01946361,
0.02071939, 0.02205619, 0.02347924, 0.02499411, 0.02660671, 0.02832336,
0.03015076, 0.03209607, 0.03416689, 0.03637131, 0.03871796, 0.04121602,
0.04387525, 0.04670605, 0.0497195, 0.05292737, 0.05634221, 0.05997737,
0.06384707, 0.06796645, 0.0723516, 0.07701967, 0.08198893, 0.08727881,
0.09290998, 0.09890447, 0.1052857, 0.1120787, 0.1193099, 0.1270077,
0.1352022, 0.1439253, 0.1532113, 0.1630964, 0.1736193, 0.1848211,
0.1967456, 0.2094395])
def display_polar_spectra(allspecs,heading,part='Re',limit_display_wv=100,title=None,outputfile=None,
fig=None,ax=None,interactive=False,tau_id=2):
for tauval in [tau_id] :
if part=='Im':
thecmap='PuOr'
coS = allspecs['cross-spectrum_%stau' % tauval].mean(
dim='%stau' % tauval).imag # co spectrum = imag part of cross-spectrum
else:
thecmap=cmap
coS = abs(allspecs['cross-spectrum_%stau' % tauval].mean(
dim='%stau' % tauval).real) # co spectrum = real part of cross-spectrum
new_spec_Polar = xsarsea.conversion_polar_cartesian.from_xCartesianSpectrum(coS,Nphi=72,
ksampling='log',**{'Nk' : 60,
'kmin' :
reference_oswK_1145m_60pts[
0],
'kmax' :
reference_oswK_1145m_60pts[
-1]})
crossSpectraRePol = new_spec_Polar.squeeze()
# crossSpectraRePol = spectrum_clockwise_to_trigo.apply_clockwise_to_trigo(crossSpectraRePol)
# crossSpectraRePol = spectrum_rotation.apply_rotation(crossSpectraRePol,
# 90.) # This is for having origin at North
crossSpectraRePol = spectrum_rotation.apply_rotation(crossSpectraRePol,heading)
# partie display plot
if fig is None:
plt.figure(figsize=(20,8),dpi=100)
if ax is None:
ax = plt.subplot(1,1,1,polar=True)
ax.set_theta_direction(-1) # theta increasing clockwise
ax.set_theta_zero_location("N")
print(np.log(abs(crossSpectraRePol)).min())
print(crossSpectraRePol.min(),crossSpectraRePol.max(),crossSpectraRePol.mean())
crossSpectraRe_redpol = crossSpectraRePol.where(np.abs(new_spec_Polar.k) <= 2 * np.pi / limit_display_wv,
drop=True)
print('heading',heading)
#shift_phi = np.radians(heading-90)
crossSpectraRe_redpol_phi_new = crossSpectraRe_redpol
#crossSpectraRe_redpol_phi_new = crossSpectraRe_redpol.assign_coords({'phi':(crossSpectraRe_redpol.phi.values+shift_phi)[::-1]})
# crossSpectraRe_redpol_phi_new = crossSpectraRe_redpol_phi_new/45000.
crossSpectraRe_redpol_phi_new.plot(cmap=thecmap,alpha=0.8)
plt.grid(True)
plt.plot(np.radians([heading,heading]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.plot(np.radians([heading + 180,heading + 180]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.plot(np.radians([heading + 90,heading + 90]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.plot(np.radians([heading + 270,heading + 270]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
ax.text(np.radians(heading),(0.65 + (np.cos(np.radians(heading)) < 0) * 0.25) * max(crossSpectraRe_redpol.k),
' Azimuth',size=18,color='red',rotation=-heading + 90 + (np.sin(np.radians(heading)) < 0) * 180,
ha='center')
ax.text(np.radians(heading + 90),0.80 * max(crossSpectraRe_redpol.k),' Range',size=18,color='red',
rotation=-heading + (np.cos(np.radians(heading)) < 0) * 180,ha='center')
display_xspectra_grid.circle_plot(ax,r=[300,400,500,600])
ax.set_rmax(2.0*np.pi/limit_display_wv)
if title is None:
plt.title('tau : %s' % tauval,fontsize=18)
else:
plt.title(title,fontsize=18)
if outputfile:
plt.savefig(outputfile)
logging.info('outputfile : %s',outputfile)
else:
if interactive:
plt.show()
def display_cartesian_spectra(allspecs,part='Re',limit_display_wv=100,title=None,outputfile=None,
fig=None,ax=None,interactive=False,tau_id=2):
for tauval in [tau_id] :
if part=='Im':
thecmap='PuOr'
coS = allspecs['cross-spectrum_%stau' % tauval].mean(
dim='%stau' % tauval).imag # co spectrum = imag part of cross-spectrum
else:
thecmap=cmap
coS = abs(allspecs['cross-spectrum_%stau' % tauval].mean(
dim='%stau' % tauval).real) # co spectrum = real part of cross-spectrum
new_spec_Polar = xsarsea.conversion_polar_cartesian.from_xCartesianSpectrum(coS,Nphi=72,
ksampling='log',**{'Nk' : 60,
'kmin' :
reference_oswK_1145m_60pts[
0],
'kmax' :
reference_oswK_1145m_60pts[
-1]})
crossSpectraRePol = new_spec_Polar.squeeze()
# crossSpectraRePol = spectrum_clockwise_to_trigo.apply_clockwise_to_trigo(crossSpectraRePol)
# crossSpectraRePol = spectrum_rotation.apply_rotation(crossSpectraRePol,
# 90.) # This is for having origin at North
crossSpectraRePol = spectrum_rotation.apply_rotation(crossSpectraRePol,heading)
# partie display plot
if fig is None:
plt.figure(figsize=(20,8),dpi=100)
if ax is None:
ax = plt.subplot(1,1,1,polar=True)
ax.set_theta_direction(-1) # theta increasing clockwise
ax.set_theta_zero_location("N")
print(np.log(abs(crossSpectraRePol)).min())
print(crossSpectraRePol.min(),crossSpectraRePol.max(),crossSpectraRePol.mean())
crossSpectraRe_redpol = crossSpectraRePol.where(np.abs(new_spec_Polar.k) <= 2 * np.pi / limit_display_wv,
drop=True)
print('heading',heading)
#shift_phi = np.radians(heading-90)
crossSpectraRe_redpol_phi_new = crossSpectraRe_redpol
#crossSpectraRe_redpol_phi_new = crossSpectraRe_redpol.assign_coords({'phi':(crossSpectraRe_redpol.phi.values+shift_phi)[::-1]})
# crossSpectraRe_redpol_phi_new = crossSpectraRe_redpol_phi_new/45000.
crossSpectraRe_redpol_phi_new.plot(cmap=thecmap,alpha=0.8)
plt.grid(True)
plt.plot(np.radians([heading,heading]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.plot(np.radians([heading + 180,heading + 180]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.plot(np.radians([heading + 90,heading + 90]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.plot(np.radians([heading + 270,heading + 270]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
ax.text(np.radians(heading),(0.65 + (np.cos(np.radians(heading)) < 0) * 0.25) * max(crossSpectraRe_redpol.k),
' Azimuth',size=18,color='red',rotation=-heading + 90 + (np.sin(np.radians(heading)) < 0) * 180,
ha='center')
ax.text(np.radians(heading + 90),0.80 * max(crossSpectraRe_redpol.k),' Range',size=18,color='red',
rotation=-heading + (np.cos(np.radians(heading)) < 0) * 180,ha='center')
display_xspectra_grid.circle_plot(ax,r=[300,400,500,600])
ax.set_rmax(2.0*np.pi/limit_display_wv)
if title is None:
plt.title('tau : %s' % tauval,fontsize=18)
else:
plt.title(title,fontsize=18)
if outputfile:
plt.savefig(outputfile)
logging.info('outputfile : %s',outputfile)
else:
if interactive:
plt.show()
......@@ -5,6 +5,8 @@ import logging
import sys
from matplotlib import colors as mcolors
from xsarsea.conversion_polar_cartesian import from_xCartesianSpectrum
import spectrum_clockwise_to_trigo
import spectrum_rotation
#sys.path.append('/home1/datahome/agrouaze/git/sar/sar')
#from utils.gmf import compute_roughness # depot sar
def draw_signal_and_grid(slc,splitting_image,one_tiff,variable_displayed='digital_number'):
......@@ -113,10 +115,11 @@ def draw_signal_and_grid(slc,splitting_image,one_tiff,variable_displayed='digita
plt.show()
def display_xspectra_per_domain(allspecs_per_sub_domain):
def display_xspectra_per_domain(allspecs_per_sub_domain,part='Re'):
"""
:param allspecs_per_sub_domain:
:param (str): Re or Im, component of complex x spectra to display
:return:
"""
cmap = mcolors.LinearSegmentedColormap.from_list("",["white","violet","mediumpurple","cyan","springgreen","yellow",
......@@ -154,19 +157,50 @@ def display_xspectra_per_domain(allspecs_per_sub_domain):
# print('subdom',subdom,'->',xindplot,yindplot)
# plt.subplot(4,4,subdom+1)
plt.title('sub domain #%s' % subdom)
crossSpectraRe = np.abs(allspecs_per_sub_domain[subdom]['cross-spectrum_2tau'].mean(dim='2tau').real)
crossSpectraRe_red = crossSpectraRe.where(
np.logical_and(np.abs(crossSpectraRe.kx) <= 0.14,np.abs(crossSpectraRe.ky) <= 0.14),drop=True)
if vmin is None :
vmin = 0
vmax = crossSpectraRe_red.max()
levels = np.linspace(vmin,vmax,25)
crossSpectraRe_red.plot(x='kx',cmap=cmap,levels=levels,vmin=0)
if part=='Re':
crossSpectraRe = np.abs(allspecs_per_sub_domain[subdom]['cross-spectrum_2tau'].mean(dim='2tau').real)
crossSpectraRe_red = crossSpectraRe.where(
np.logical_and(np.abs(crossSpectraRe.kx) <= 0.14,np.abs(crossSpectraRe.ky) <= 0.14),drop=True)
crossSpectraRe_red = crossSpectraRe_red.rolling(kx=3,center=True).mean().rolling(ky=3,center=True).mean()
if vmin is None :
vmin = 0
vmax = crossSpectraRe_red.max()
levels = np.linspace(vmin,vmax,25)
crossSpectraRe_red.plot(x='kx',cmap=cmap,levels=levels,vmin=0)
else:
crossSpectraIm = allspecs_per_sub_domain[subdom]['cross-spectrum_2tau'].mean(dim='2tau').imag
crossSpectraIm_red = crossSpectraIm.where(
np.logical_and(np.abs(crossSpectraIm.kx) <= 0.14,np.abs(crossSpectraIm.ky) <= 0.14),drop=True)
crossSpectraIm_red = crossSpectraIm_red.rolling(kx=3,center=True).mean().rolling(ky=3,center=True).mean()
if vmin is None :
vmin = -crossSpectraIm_red.max()
vmin = crossSpectraIm_red.min()
vmax = crossSpectraIm_red.max()
levels = np.linspace(vmin,vmax,25)
crossSpectraIm_red.plot(x='kx',cmap=PuOr,levels=levels,vmin=vmin)
plt.grid()
plt.axis('scaled')
def display_xspectra_per_domain_polar(allspecs_per_sub_domain,heading):
def circle_plot ( ax,r,freq=0 ) :
# Radial Circles and their label
theta = 2 * np.pi * np.arange(360) / 360.
labels = []
if freq == 0 :
for i in r :
plt.plot(theta,np.arange(360) * 0 + 2 * np.pi / i,'--k')
labels.append(str(i) + ' m')
ax.set_rgrids([2 * np.pi / i for i in r],labels=labels,angle=45.)
if freq == 1 :
for i in r :
plt.plot(theta,np.arange(360) * 0 + np.sqrt(2 * np.pi / i * 9.81 / (2 * np.pi) ** 2),'--k')
labels.append(str(i) + ' m')
ax.set_rgrids([np.sqrt(2 * np.pi / i * 9.81 / (2 * np.pi) ** 2) for i in r],labels=labels,angle=45.)
def display_xspectra_per_domain_polar(allspecs_per_sub_domain,heading,part='Re',min_wavelength=100):
"""
:param allspecs_per_sub_domain:
......@@ -200,15 +234,25 @@ def display_xspectra_per_domain_polar(allspecs_per_sub_domain,heading):
xindplot,yindplot = np.where(final_inds == subdom)
# print()
ax = fig.add_subplot(specgrid[xindplot[0],yindplot[0]],polar=True)
ax.set_theta_direction(-1) # theta increasing clockwise
ax.set_theta_zero_location("N")
# print('subdom',subdom,'->',xindplot,yindplot)
# plt.subplot(4,4,subdom+1)
plt.title('sub domain #%s' % subdom)
crossSpectraRe = np.abs(allspecs_per_sub_domain[subdom]['cross-spectrum_2tau'].mean(dim='2tau').real)
if part=='Re':
crossSpectraRe = np.abs(allspecs_per_sub_domain[subdom]['cross-spectrum_2tau'].mean(dim='2tau').real)
else:
crossSpectraRe = allspecs_per_sub_domain[subdom]['cross-spectrum_2tau'].mean(dim='2tau').imag
#crossSpectraRe = crossSpectraRe.rolling(kx=3,center=True).mean().rolling(ky=3,center=True).mean()
logging.warning('trick to match a reference spectra on cartesian coords kx and ky')
cspRe_okGrid = crossSpectraRe
#cspRe_okGrid = crossSpectraRe
# cspRe_okGrid = crossSpectraRe.assign_coords(kx=crossSpectraRe.kx.data / (1.1 * np.pi),
# ky=crossSpectraRe.ky.data / (1.1 * np.pi))
new_spec_PolarRe = from_xCartesianSpectrum(cspRe_okGrid,Nphi=72,ksampling='log',**{'Nk' : 60})
crossSpectraRe = from_xCartesianSpectrum(crossSpectraRe,Nphi=72,ksampling='log',**{'Nk' : 60})
crossSpectraRe = spectrum_clockwise_to_trigo.apply_clockwise_to_trigo(crossSpectraRe)
crossSpectraRe = spectrum_rotation.apply_rotation(crossSpectraRe,
90.) # This is for having origin at North
crossSpectraRe = spectrum_rotation.apply_rotation(crossSpectraRe,heading)
# crossSpectraRe_red = crossSpectraRe.where(
# np.logical_and(np.abs(crossSpectraRe.kx) <= 0.14,np.abs(crossSpectraRe.ky) <= 0.14),drop=True)
# ax.set_rmax(0.8)
......@@ -219,22 +263,31 @@ def display_xspectra_per_domain_polar(allspecs_per_sub_domain,heading):
#print(crossSpectraRePol.min(),crossSpectraRePol.max(),crossSpectraRePol.mean())
# plt.pcolor(crossSpectraRePol/crossSpectraRePol.max())
# plt.contourf(new_spec_Polar.phi,new_spec_Polar.k,crossSpectraRePol)
crossSpectraRe_redpol = new_spec_PolarRe.where(np.abs(new_spec_PolarRe.k) <= 0.05,drop=True)
if vmin is None :
vmin = 0
crossSpectraRe_redpol = crossSpectraRe.where(np.abs(crossSpectraRe.k) <= 2.*np.pi/min_wavelength,drop=True)
crossSpectraRe_redpol = crossSpectraRe_redpol.rolling(k=3,center=True).mean().rolling(phi=3,center=True).mean()
if part=='Re':
if vmin is None :
vmin = 0
else:
vmin = crossSpectraRe_redpol.min()
vmax = crossSpectraRe_redpol.max()
levels = np.linspace(vmin,vmax,25)
crossSpectraRe_redpol.plot(cmap=cmap,levels=levels,vmin=0)
plt.polar(np.radians([heading,heading]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.polar(np.radians([heading + 180,heading + 180]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.polar(np.radians([heading + 90,heading + 90]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
plt.polar(np.radians([heading + 270,heading + 270]),(0.01,max(crossSpectraRe_redpol.k)),'r-')
ax.text(np.radians(heading),(0.65 + (np.cos(np.radians(heading)) < 0) * 0.25) * max(crossSpectraRe_redpol.k),
if part == 'Re' :
crossSpectraRe_redpol.plot(cmap=cmap,levels=levels,vmin=0)
else:
crossSpectraRe_redpol.plot(cmap=PuOr,levels=levels)
max_k = max(crossSpectraRe_redpol.k.values)
#print('max_k',max_k)
plt.plot(np.radians([heading,heading]),(0.01,max_k),'r-')
plt.plot(np.radians([heading + 180,heading + 180]),(0.01,max_k),'r-')
plt.plot(np.radians([heading + 90,heading + 90]),(0.01,max_k),'r-')
plt.plot(np.radians([heading + 270,heading + 270]),(0.01,max_k),'r-')
ax.text(np.radians(heading),(0.65 + (np.cos(np.radians(heading)) < 0) * 0.25) * max_k,
' Azimuth',size=18,color='red',rotation=-heading + 90 + (np.sin(np.radians(heading)) < 0) * 180,
ha='center')
ax.text(np.radians(heading + 90),0.80 * max(crossSpectraRe_redpol.k),' Range',size=18,color='red',
ax.text(np.radians(heading + 90),0.80 * max_k,' Range',size=18,color='red',
rotation=-heading + (np.cos(np.radians(heading)) < 0) * 180,ha='center')
plt.grid(0,axis='y')
plt.axis('scaled')
\ No newline at end of file
r = [200,300,400,600,1000]
circle_plot(ax,r,freq=0)
#plt.axis('scaled')
\ No newline at end of file
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