Commit 8b364b3e authored by mickael's avatar mickael

add buoy track from chuan

parent b1728bca
255 255 255
255 255 255
241 233 250
228 211 245
214 190 240
201 168 235
187 146 230
174 125 225
161 109 220
149 93 215
137 78 210
125 62 205
113 46 200
101 31 195
89 15 190
77 0 186
72 3 190
67 6 194
62 10 198
57 13 203
52 17 207
48 20 211
43 24 216
38 27 220
33 30 224
28 34 229
24 37 233
19 41 237
14 44 242
9 48 246
4 51 250
0 55 255
0 62 255
0 70 255
0 78 255
0 86 255
0 94 255
0 102 255
0 110 255
0 118 255
0 126 255
0 134 255
0 142 255
0 150 255
0 158 255
0 166 255
0 174 255
0 182 255
0 185 255
0 189 255
0 193 255
0 197 255
0 200 255
0 204 255
0 208 255
0 212 255
0 215 255
0 219 255
0 223 255
0 227 255
0 230 255
0 234 255
0 238 255
0 242 255
0 242 251
0 242 248
0 243 245
0 243 242
0 244 238
0 244 235
0 244 232
0 245 229
0 245 226
0 246 222
0 246 219
0 246 216
0 247 213
0 247 209
0 248 206
0 248 203
0 248 200
0 249 197
0 249 193
0 250 190
0 250 187
0 250 184
0 251 180
0 251 177
0 252 174
0 252 171
0 252 168
0 253 164
0 253 161
0 254 158
0 254 155
0 255 152
4 255 147
8 255 142
12 255 137
16 255 133
20 255 128
24 255 123
28 255 118
32 255 114
36 255 109
40 255 104
45 255 99
49 255 95
53 255 90
57 255 85
61 255 80
65 255 76
69 255 71
73 255 66
77 255 61
81 255 57
85 255 52
90 255 47
94 255 42
98 255 38
102 255 33
106 255 28
110 255 23
114 255 19
118 255 14
122 255 9
126 255 4
131 255 0
134 255 0
138 255 0
142 255 0
146 255 0
150 255 0
154 255 0
158 255 0
162 255 0
165 255 0
169 255 0
173 255 0
177 255 0
181 255 0
185 255 0
189 255 0
193 255 0
196 255 0
200 255 0
204 255 0
208 255 0
212 255 0
216 255 0
220 255 0
224 255 0
227 255 0
231 255 0
235 255 0
239 255 0
243 255 0
247 255 0
251 255 0
255 255 0
255 252 0
255 249 0
255 246 0
255 243 0
255 240 0
255 238 0
255 235 0
255 232 0
255 229 0
255 226 0
255 224 0
255 221 0
255 218 0
255 215 0
255 212 0
255 210 0
255 207 0
255 204 0
255 201 0
255 198 0
255 195 0
255 193 0
255 190 0
255 187 0
255 184 0
255 181 0
255 179 0
255 176 0
255 173 0
255 170 0
255 167 0
255 165 0
255 162 0
255 160 0
255 158 0
255 156 0
255 153 0
255 151 0
255 149 0
255 147 0
255 145 0
255 142 0
255 140 0
255 138 0
255 136 0
255 133 0
255 131 0
255 129 0
255 127 0
255 125 0
255 122 0
255 120 0
255 118 0
255 116 0
255 113 0
255 111 0
255 109 0
255 107 0
255 105 0
255 102 0
255 100 0
255 98 0
255 96 0
255 94 0
255 88 0
255 83 0
255 77 0
255 72 0
255 66 0
255 61 0
255 55 0
255 50 0
255 45 0
255 39 0
255 34 0
255 28 0
255 23 0
255 17 0
255 12 0
255 7 0
251 6 0
248 6 0
245 5 0
242 5 0
239 4 0
236 4 0
233 3 0
230 3 0
227 2 0
224 2 0
221 1 0
218 1 0
215 0 0
212 0 0
209 0 0
0 0 0
This diff is collapsed.
This diff is collapsed.
import cPickle,bz2
def zip_save(filename,*objects):
with bz2.BZ2File(filename, 'wb', compresslevel=9) as f:
# with gzip.GzipFile(filename, 'wb') as f:
for obj in objects: cPickle.dump(obj, f,2)
def zip_load(filename):
with bz2.BZ2File(filename, 'rb', compresslevel=9) as f:
# with gzip.GzipFile(filename, 'rb') as f:
data = cPickle.load(f)
return data
\ No newline at end of file
ndbc_globwave_Path = '/home/cercache/project/globwave/data/globwave/ndbc/archive/wave_sensor/'
cdip_globwave_Path = '/home/cercache/project/globwave/data/globwave/cdip/wave_sensor/'
ndbc_nrt_Path = '/home/cercache/project/globwave/data/provider/ndbc/spectra/nrt/'
archive_globwave = '/home/cercache/project/globwave/data/globwave/ndbc/archive/wave_sensor/'
ndbc_archive_Path = '/home/cercache/project/globwave/data/provider/ndbc/spectra/'
cdip_txt_Path = '/home/cercache/project/globwave/data/provider/cdip/'
# -*- coding: utf-8 -*-
# pour decoder la date, voir float2strdate
import netCDF4 as nc
from glob import glob
import numpy as np
from datetime import *
import pdb
def get_oswl2ocn_param(file,get_spec=False,norut=None):
data = nc.Dataset(file)
#pdb.set_trace()
wv_param={}
wv_param['oswHs'] = (data.variables['oswHs'][:])
wv_param['oswWl'] = (data.variables['oswWl'][:])
wv_param['oswWindSeaHs'] = (data.variables['oswWindSeaHs'][:])
wv_param['oswAzCutoff'] = (data.variables['oswAzCutoff'][:])
wv_param['oswRaCutoff'] = (data.variables['oswRaCutoff'][:])
wv_param['oswSnr'] = (data.variables['oswSnr'][:])
wv_param['oswNv'] = (data.variables['oswNv'][:])
wv_param['oswInten'] = (data.variables['oswInten'][:])
wv_param['oswWindSpeed'] = (data.variables['oswWindSpeed'][:])
wv_param['oswWindDirection'] = (data.variables['oswWindDirection'][:])
wv_param['oswIncidenceAngle'] = (data.variables['oswIncidenceAngle'][:])
wv_param['oswLandFlag'] = (data.variables['oswLandFlag'][:])
wv_param['oswDepth'] = (data.variables['oswDepth'][:])
wv_param['oswLandCoverage'] = (data.variables['oswLandCoverage'][:])
wv_param['oswLat'] = (data.variables['oswLat'][:])
wv_param['oswLon'] = (data.variables['oswLon'][:])
wv_param['oswIconf'] = (data.variables['oswIconf'][:])
wv_param['oswDirmet'] = (data.variables['oswDirmet'][:])
wv_param['oswNrcs'] = (data.variables['oswNrcs'][:])
wv_param['oswHeading'] = (data.variables['oswHeading'][:])
wv_param['oswSpecRes'] = (data.variables['oswSpecRes'][:][0,0,:])
wv_param['oswWaveAge'] = (data.variables['oswWaveAge'][:])
#wv_param['rvlDcObs'] = (data.variables['rvlDcObs'][:])
#wv_param['rvlDcGeo'] = (data.variables['rvlDcGeo'][:])
#wv_param['rvlDcMiss'] = (data.variables['rvlDcMiss'][:])
#try:
# wv_param['oswEcmwfWindDirection'] = (data.variables['oswEcmwfWindDirection'][:])
#except KeyError:
# wv_param['oswEcmwfWindDirection'] = np.nan
if get_spec==True:
if norut is None:
wv_param['oswk'] = (data.variables['oswK'][:])
wv_param['oswphi'] = (data.variables['oswPhi'][:])
wv_param['oswspec'] = (data.variables['oswPolSpec'][0,0,:,:])
wv_param['xspec_re'] = (data.variables['oswQualityCrossSpectraRe'][0,0,:,:])
wv_param['xspec_im'] = (data.variables['oswQualityCrossSpectraIm'][0,0,:,:])
wv_param['oswPartitions'] = (data.variables['oswPartitions'][0,0,:,:])
else:
wv_param['oswk'] = (data.variables['oswK'][:])
wv_param['oswphi'] = (data.variables['oswPhi'][:])
wv_param['oswspec'] = (data.variables['oswPolSpec'][:,0,:,:])
wv_param['xspec_re'] = (data.variables['oswQualityCrossSpectraRe'][:,0,:,:])
wv_param['xspec_im'] = (data.variables['oswQualityCrossSpectraIm'][:,0,:,:])
wv_param['oswPartitions'] = (data.variables['oswPartitions'][:,0,:,:])
#pdb.set_trace()
data.close()
#pdb.set_trace()
return wv_param
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 26 10:55:29 2017
@author: xwang
Update by Mlebars
Exemple of command to launch
python main2.py --provider ndbc_nrt --ids 32012,46041 --startdate 20170101-080000 --enddate 20170102-080000
"""
#python main2.py --provider ndbc_nrt --ids 32012,46041 --startdate 20170101-080000 --enddate 20170102-080000
import matplotlib
import numpy as np
from bz_pick import zip_load,zip_save
import matplotlib as mpl
import matplotlib.pyplot as plt
from datetime import datetime,timedelta
from buoy_nc_reader import read_buoy_nc,read_nc_wind,spectrum_mean
import pdb
from data_path import cdip_globwave_Path, ndbc_globwave_Path,ndbc_nrt_Path,ndbc_archive_Path
from palette import *
from buoy_and_SAR_storm_plot import buoy_partition_plot, plot_Partition
from glob import glob
#MLB Correction
import os
import argparse
import sys
Method='G'
g=9.8
mcolor=['b','m','r','w','g','c','y']
#provider='ndbc_globwave'
provider = 'ndbc_nrt'
dirfig='/home/cercache/users/mlebars/buoyspec/'+provider+'/pic_swell_G/'
dir_output ='/home/cercache/users/mlebars/buoyspec/'+provider+'/BuoySwell_G/'
def cmdlineparser():
"""
define the options of the scripts and set func depending of the inputs
"""
#OBJCHOICES = ('WV', 'SAFE')
parser = argparse.ArgumentParser(description='Sentinel L2 Indexer')
parser.add_argument('--provider', action='store', default=None,choices=['ndbc_nrt','ndbc_globwave'],required=True, help='Mode that determine to search data')
parser.add_argument('--ids', action='store',type=str, default=None, help='give the buoys ids ex: 32012,46041,46076',required=True)
parser.add_argument('--startdate', action='store',type=str, default=None, help='startDate to process YYYYMMDD-hhmmss',required=True)
parser.add_argument('--enddate', action='store',type=str, default=None, help='endDate to process YYYYMMDD-hhmmss',required=True)
#parser.add_argument('data_type', type=str.upper, choices=OBJCHOICES)
#subparsers = parser.add_subparsers()
#parser_wv = subparsers.add_parser('wv')
#parser_wv.set_defaults(func=run_l2_wv_indexation)
#parser_wv_phx = subparsers.add_parser('wvphx')
#parser_wv_phx.set_defaults(func=run_l2_wv_indexation_phoenix)
return parser
def main():
parser = cmdlineparser()
args = parser.parse_args()
ids = args.ids.split(',')
time_start=datetime.strptime(args.startdate, '%Y%m%d-%H%M%S')
time_end=datetime.strptime(args.enddate, '%Y%m%d-%H%M%S')
if Method=='G':
print time_start,time_end
print "started at",datetime.now()
#MODIF MLB
#from waveCluster_buoy_G import track_buoy,plot_buoy_partitions_cluster,create_buoy_nc
from waveCluster_buoy_G_try import track_buoy,plot_buoy_partitions_cluster,create_buoy_nc
for id in ids:
print 'Buoy Id:', id
time_wd = {'type':'t1t2','startDate': time_start-timedelta(days = 1),'stopDate':time_end+timedelta(days = 1)}
phi = np.arange(0,360,10)
lon,lat,time,f,k,sf,sp_list,buoy_flag = read_buoy_nc(id,time_wd,phi = phi,con = 'to', provider = provider,smooth=True)
if lon is None:
print "lon not Found"
continue
wtime,wspd,wdir = read_nc_wind(id,time_wd)
_lon = np.mean(lon) % 360
_lat = np.mean(lat)
#dirfig='//home/cercache/users/xwang/buoyspec/pic_swell_G/'
if os.path.isdir(dirfig)==False:
command_line = 'mkdir ' + dirfig
os.system(command_line)
mytitle=id+' lon: {:3.2f}'.format(_lon)+' lat: {:3.2f}'.format(_lat)
fileout = dirfig + id +'_'+datetime.strftime(time_start,'%Y-%m-%d %H:%M:%S')+'buoy_partitions.png'
#Modif MLB
#[hs_buoy,wl_buoy,dir_buoy,time_buoy]=buoy_partition_plot(time,f,k,sf,sp_list,mytitle,fileout,buoy_flag,myTick='day')
[hs_buoy,wl_buoy,dir_buoy,time_buoy,parts_buoy]=buoy_partition_plot(time,f,k,sf,sp_list,mytitle,fileout,buoy_flag,myTick='day')
print "after buoy_partition_plot",fileout
if len(time_buoy)==0:
continue
#Modif MLB
#[btimelist, bfplist, bhslist, bdirlist, bsyslist,bwindlist,thr]=track_buoy(hs_buoy,wl_buoy,dir_buoy,time_buoy,id,time,f,k,sf,sp_list,wtime,wspd,wdir)
[btimelist, bfplist, bhslist, bdirlist, bsyslist,bwindlist,thr,partlist]=track_buoy(hs_buoy,wl_buoy,dir_buoy,time_buoy,parts_buoy,id,time,f,k,sf,sp_list,wtime,wspd,wdir,Method)
if max(bsyslist)<1:
continue
plot_buoy_partitions_cluster( btimelist, bfplist, bhslist, bdirlist, bsyslist,bwindlist,time,f,k,sf,sp_list,mcolor,dirfig + id,mytitle,wtime,wspd,wdir)
#dir_output ='//home/cercache/users/xwang/buoyspec/BuoySwell_G/'
#Modif MLB
#create_buoy_nc( btimelist, bfplist, bhslist, bdirlist, bsyslist,bwindlist,thr,wtime,wspd,wdir,id,_lon,_lat,dir_output)
create_buoy_nc( btimelist, bfplist, bhslist, bdirlist, bsyslist,bwindlist,thr,partlist,wtime,wspd,wdir,id,_lon,_lat,dir_output)
print (dir_output)
print datetime.now()
print "End",datetime.now()
sys.exit(0)
if __name__ == '__main__':
main()
def myFunction():
#archivePath = ndbc_globwave_Path
#filepattern = 'WMO*spectrum*.nc'
#monthDir = os.path.join(archivePath, str(time_end.year), "%02d" % time_end.month)
#files = glob(monthDir+'/'+filepattern)
ids = []
#for _fi in files:
# id = str.split(str.split(os.path.basename(_fi),'_')[0],'WMO')[1]
# ids.append(id)
ids = ['32012','46041','46076','46002','46042','46078','46006','46047','46083','46011','46053','46084','46012','46054','46086','46014','46059','46087','46015','46060','46089','46022','46061','51000','46026','46066','51003','46027','46069','51004','46028','46070','51101','46029','46071','46035','46072']
#ids = ['32012']
\ No newline at end of file
This diff is collapsed.
import numpy as np
from numpy import sin,cos,pi
import pdb
def buoy_spectrum2d(sf,a1,a2,b1,b2,f,smooth,dirs = np.arange(0,365,10)):
# Maximum entropy method to estimate the Directional Distribution
# Maximum Entropy Method - Lygre & Krogstad (1986 - JPO)
# Eqn. 13:
# phi1 = (c1 - c2c1*)/(1 - abs(c1)^2)
# phi2 = c2 - c1phi1
# 2piD = (1 - phi1c1* - phi2c2*)/abs(1 - phi1exp(-itheta) -phi2exp(2itheta))^2
# c1 and c2 are the complex fourier coefficients
# inputs : dirs in degrees
nfreq = sf.size
nbin = dirs.size
c1 = a1+1j*b1
c2 = a2+1j*b2
p1 = (c1-c2*c1.conjugate())/(1.-abs(c1)**2)
p2 = c2-c1*p1
# numerator(2D) : x
x = 1.-p1*c1.conjugate()-p2*c2.conjugate()
x= np.tile(np.real(x),(nbin,1)).T
# denominator(2D): y
a = dirs*pi/180.
e1 = np.tile(cos(a)-1j*sin(a),(nfreq,1))
e2 = np.tile(cos(2*a)-1j*sin(2*a),(nfreq,1))
p1e1 = np.tile(p1,(nbin,1)).T*e1
p2e2 = np.tile(p2,(nbin,1)).T*e2
y = abs(1-p1e1-p2e2)**2
D = x/(y*2*pi)
# normalizes the spreading function,
# so that int D(theta,f) dtheta = 1 for each f
dth = dirs[1]-dirs[0]
tot = np.tile(np.sum(D, axis=1)*dth/180.*pi,(nbin,1)).T
D = D/tot
sp2d = np.tile(sf,(nbin,1)).T*D
# sp2d=np.nan_to_num(sp2d)
# smooth
if smooth==True:
from scipy.signal import fftconvolve
df=[]
df.append(f[1]-f[0])
for ind in range(len(f)-2):
df.append(1.0/2*(f[ind+2]-f[ind]))
df.append(f[-1]-f[-2])
df_grid = np.tile(df, (len(dirs), 1))
# pdb.set_trace()
sp2d=sp2d*df_grid.T
K = np.array([[1./np.sqrt(2),1,1./np.sqrt(2)],[1,2,1],[1./np.sqrt(2),1,1./np.sqrt(2)]])/(6+4./np.sqrt(2))
sp2d = fftconvolve(sp2d,K,mode = 'same')
sp2d = sp2d*(1.0/df_grid.T)
# pdb.set_trace()
return sp2d,D
def spfphi2kphi(spfphi,f):
g = 9.81
k = (2*pi*f)**2/g
dk = np.zeros(k.size)
dk[0] = k[1]-k[0]
dk[-1] = k[-1]-k[-2]
dk[1:-1] = (k[2:]-k[0:-2])/2.
df = np.zeros(f.size)
df[0] = f[1]-f[0]
df[-1] = f[-1]-f[-2]
df[1:-1] = (f[2:]-f[0:-2])/2.
spkphi = []
for sp in spfphi:
df2kdk = (np.tile(df/(k*dk), (sp.shape[1], 1))).T
sp2 = sp*df2kdk
spkphi.append(sp2)
return spkphi,k
def spfrom2to(sp_in,phi):
phi2 = (180.+phi)%360
ind = np.argsort(phi2)
sp_out = []
for sp in sp_in:
sp_out.append(sp[:,ind])
return sp_out
def th2ab(stheta1,stheta2,theta1,theta2,construct=None):
# convert theta & stheta(globwave data used) to a1,b1,a2,b2
r1 = abs(1-0.5*(stheta1/180.*pi)**2)
r2 = abs(1-2.0*(stheta2/180.*pi)**2)
# r1 = 1-0.5*(stheta1/180.*pi)**2
# r2 = 1-2.0*(stheta2/180.*pi)**2
a1 = r1*cos(theta1/180.*pi)
b1 = r1*sin(theta1/180.*pi)
a2 = r2*cos(2*theta2/180.*pi)
b2 = r2*sin(2*theta2/180.*pi)
return a1,a2,b1,b2
def ar2ab(alpha1,alpha2,r1,r2):
# convert alpha & r(ndbc data used) to a1,b1,a2,b2
a1 = r1*cos(alpha1/180.*pi)
b1 = r1*sin(alpha1/180.*pi)
a2 = r2*cos(2*alpha2/180.*pi)
b2 = r2*sin(2*alpha2/180.*pi)
return a1,a2,b1,b2
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
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