forked from rdkreij/oceanrow
-
Notifications
You must be signed in to change notification settings - Fork 0
/
mod_ocean_row.py
1410 lines (1199 loc) · 53.8 KB
/
mod_ocean_row.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# import packages
import pandas as pd
import numpy as np
import cartopy.io.shapereader as shpreader
import matplotlib as mpl
import matplotlib.pyplot as plt
import xarray as xr
from geographiclib.geodesic import Geodesic
from shapely import geometry
from shapely.ops import unary_union
from descartes import PolygonPatch
from scipy.optimize import fsolve,root
from shapely.geometry import Point
import calendar
import cartopy.crs as ccrs
import cartopy.feature as cfeature
class ocean_row:
'''
Class to simulate ocean rows
'''
def __init__(self,mcf):
'''
Initialize
self
----
mcf: dictionary
start_lat: scalar
starting location latitude [deg]
start_lon: scalar
starting location longitude [deg]
start_name: scalar
starting location name
stop_lat: scalar
final destination latitude [deg]
stop_lon: scalar
final destination longitude [deg]
stop_name: scalar
final destination name
vpers: scalar
velocity in perfect conditions [m/s]
tod_0: datetime
starting timestamp
tod_start: list of strings ('HH:MM:SS')
list of timestamps to start rowing
tod_stop: list of strings ('HH:MM:SS')
list of timestamps to stop rowing
m: scalar
mass
Gw: scalar
water friction constant
GwA: scalar
water friction constant sea anchor
Ga: scalar
air friction constant
'''
self.mcf = mcf
def distance(self):
'''
Calculate distance between starting and final destination
Parameters
----------
self
Returns
-------
distance: scalar
distance between points [m]
'''
lat1 = self.mcf['start_lat']
lon1 = self.mcf['start_lon']
lat2 = self.mcf['stop_lat']
lon2 = self.mcf['stop_lon']
return Geodesic.WGS84.Inverse(lat1,lon1,lat2,lon2)['s12']
def rowing_force(self):
'''
Calculate rowing force given the velocity in perfect conditions
Parameters
----------
self
Returns
-------
Frow: scalar
rowing force
'''
Ga = self.mcf['Ga']
Gw = self.mcf['Gw']
vperf = self.mcf['vperf']
return (Ga+Gw)*vperf**2
def rowing_activity(self):
'''
Construct intervals indicating rowing activity
Parameters
----------
self
Returns
-------
int0:
interval starting times [s of day]
introw:
indicate per interval whether rowing (True) or not rowing (False) [True/False]
intdur:
duration per interval [s]
intN: total number of intervals per day
k0: interval index
dt0: timestep till next interval
'''
tod_0 = self.mcf['tod_0']
tod_start = self.mcf['tod_start']
tod_stop = self.mcf['tod_stop']
int0row = np.array(pd.to_timedelta(tod_start).total_seconds()).astype(int) # start rowing [s]
int0rest = np.array(pd.to_timedelta(tod_stop).total_seconds()).astype(int) # start resting [s]
int0 = np.append(int0row,int0rest) # interval starting times (both rowing and resting) [s]
introw = np.append(np.ones(int0row.shape).astype(bool),np.zeros(int0rest.shape).astype(bool)) # activity (True = rowing, False = resting)
order = np.argsort(int0) # sort intervals in time
int0 = int0[order] # sort starting times [s]
introw = introw[order] # sort activity [bool]
intdur = np.diff(np.append(int0,24*3600+int0[0])) # calculate interval duration [s]
intN = len(int0) # total number of intervals per day
sod_t0 = int(pd.to_timedelta(tod_0[11:19]).total_seconds()) # starting time [second of day]
if sod_t0<int0[0]: # if in interval that started the day before and runs in the next day
k0 = intN-1 # interval index
dt0 = int0[0]-sod_t0 # timestep till next interval
else:
k0 = np.where(int0<=sod_t0)[0][-1] # interval index
dt0 = intdur[k0]-(sod_t0-int0[k0]) # timestep till next interval
return int0,introw,intdur,intN,k0,dt0
def plt_rowing_activity(self):
'''
Plot rowing activity as a function of time
Parameters
----------
self:
tod_0: starting timestamp [datetime]
tod_start: list of timestamps to start rowing
tod_stop: list of timestamps to stop rowing
Returns
-------
fig: created figure
ax: axes of created figure
'''
int0,introw,intdur,intN,k0,dt0 = self.rowing_activity() # row activity
# plot rowing activity as function of time of day
fig,ax = plt.subplots()
plt.step(np.append(np.append(0,int0),3600*24)/(3600),np.append(np.append(introw[-1],introw),introw[-1]),where='post',color='black')
ax.set_xlim([0,24])
ax.set_ylabel('Rowing activity')
ax.set_xlabel('$t$ (h)')
plt.show()
plt.close()
return fig,ax
@staticmethod
def plt_modified_imshow(X,Y,Z,ax=None,cb=False,*args,**kwargs): # modified plt.imshow
'''
Modified plt.imshow that automatically includes colorbar
Parameters
----------
X,Y: x[Nx],y[Ny] OR meshgrid of x and y values [Nx*Ny]
Z: data of interest 2D matrx [Nx*Ny]
ax: ploting axis
cb: include colorbar when true
args: arguments for plt.imshow
kwargs: kwargs for plt.imshow
Returns
-------
im: modified plt.imshow
'''
if X.ndim == 2:
X = X[0,:]
Y = Y[:,0]
# draw modified imshow
dx = (X[1]-X[0])/2 # spacing along x
dy = (Y[1]-Y[0])/2 # spacing along y
extent = [X[0]-dx, X[-1]+dx, Y[0]-dy, Y[-1]+dy] # extent imshow to correct size
if ax is None: # if no ax is included
ax = plt.gca() # create axis
im = ax.imshow(np.flip(Z,axis=0), extent=extent,*args,**kwargs) # draw imshow
# plot colorbar
if cb:
# data min and max
zmin = np.nanmin(Z) # minimum of data
zmax = np.nanmax(Z) # maximum of data
# get vmin value
if 'vmin' in kwargs.keys(): # check if vmin is stated
vmin = kwargs['vmin'] # get vmin value
else: # if not stated
vmin = zmin # set vmin from data
# get vmax value
if 'vmax' in kwargs.keys(): # check if vmax is stated
vmax = kwargs['vmax'] # get vmax value
else: # if not stated
vmax = zmax # set vmax from data
# determine extension of colorbar
if (zmin>=vmin) & (zmax<=vmax):
extend = 'neither'
elif (zmax>vmax) & (zmin<vmin):
extend = 'both'
elif zmax>vmax:
extend = 'max'
elif zmin<vmin:
extend = 'min'
# create colorbar on new axis
plt.colorbar(im,extend=extend,fraction = 0.047*len(Y)/len(X),pad=0.04)
return im
def plt_overview(self):
'''
Create overview plot of rowing domain and shortest distance between starting point and destination
Parameters
----------
self
Returns:
fig: fig
created figure
ax: ax
axes of created figure
'''
start_lon = self.mcf['start_lon']
start_lat = self.mcf['start_lat']
start_name = self.mcf['start_name']
stop_lon = self.mcf['stop_lon']
stop_lat = self.mcf['stop_lat']
stop_name = self.mcf['stop_name']
pdom = ocean_row.create_domain() # create allowed domain to row
# create figure of domain, including start and ending destinations
xlims = [30,130] # longitude limits [deg]
ylims = [-40,0] # latitude limits [deg]
fig,ax = ocean_row.plt_base_simple(xlims,ylims,10,10) # create figure
[ax.add_patch(PolygonPatch(pdom.geoms[0], alpha=0.2)) for i in range(len(pdom.geoms))] # plot rowing domain
# plot starting and destination locations
plt.scatter(start_lon,start_lat,color='blue') # add starting locations
plt.scatter(stop_lon,stop_lat,color='red') # add ending locations
plt.plot([start_lon,stop_lon],[start_lat,stop_lat],color='grey',zorder=0,transform=ccrs.Geodetic()) # plot line
ax.set(title=f'From {start_name} to {stop_name}, {int(self.distance()/1000)} km')
plt.show()
return fig,ax
@staticmethod
def create_domain():
'''
Create allowed domain to row
Returns
-------
pdom: Shapely MultiPolygon
Allowed rowing domain
'''
# define domain
dom = np.array([[40 ,-38 ],\
[35 ,-30 ],\
[35 ,-2 ],\
[95 ,-2 ],\
[102 ,-8 ],\
[113 ,-10 ],\
[122.5,-13.5],\
[122.5,-30 ],\
[115 ,-34 ],\
[110 ,-38 ]]) # coordinates domain (lon,lat) [deg]
pdom = geometry.Polygon([[p[0], p[1]] for p in dom]) # create polygon
# remove certain zones
dom_excl = np.array([[49.6 ,-15.2],\
[50.42,-16.1],\
[49.84,-17.2],\
[49.11,-17.19]]) # small cove Madagascar
pdom_excl = geometry.Polygon([[p[0], p[1]] for p in dom_excl]) # create polygon
pdom = pdom.difference(pdom_excl) # remove from domain
# remove land from zones
land_shp_fname = shpreader.natural_earth(resolution='50m',category='physical', name='land') # load in land
land_geom = unary_union(list(shpreader.Reader(land_shp_fname).geometries())) # get geometry
pdom = pdom.difference(land_geom) # remove from domain
return pdom
@staticmethod
def total_force(V,*arg):
'''
Calculate sum of forces on rowing boat
Parameters
----------
V: vector [2]
boat velocity
arg: list of arguments:
Vw: velocity water [longitudinal,lateral] [m/s]
Va: velocity air [longitudinal,lateral] [m/s]
ed: boat orientation [longitudinal,lateral] [-]
Frmag: rowing force [N]
row: True when rowing [True/False]
anc: Ture when anchor is dropped[True/False]
Gw: friction constant of boat with water
GwA: friction constant of anchor with water
Ga: friction constant of boat with air
Returns
-------
Ft: scalar
total force [longitudinal,lateral] [N]
'''
Vw,Va,ed,Frmag,row,anc,Gw,GwA,Ga = arg # fill in arguments
# rowing force
Fr = int(row)*Frmag*ed
# water friction sea anchor
Fwterm = np.sqrt((Vw-V).dot(Vw-V))*(Vw-V)
Fw = Gw*Fwterm # water friction
Fwa = int(anc)*GwA*Fwterm # sea anchor water friction, anchor is in when resting (row = 0)
# air friction
Fa = Ga*np.sqrt((Va-V).dot(Va-V))*(Va-V)
Fap = Fa.dot(ed)*ed # parallel component
# total force
Ft = Fw + Fwa + Fap + Fr
return Ft
@staticmethod
def tick_placer(ax,xlims,ylims,xtickspacing,ytickspacing):
'''
Place ticks in cartopy plot
Parameters
----------
ax: ax
axis of plot
xlims: list [2]
x-axis limits [deg]
xlims: list [2]
y-axis limits [deg]
xtickspacing: scalar
x tick spacing [deg]
ytickspacing: scalar
y tick spacing [deg]
Returns
-------
ax: ax
axis including placed limits
'''
# place x-ticks
xticks,xticklabels = ocean_row.tick_single_ax(xlims,xtickspacing) # tick locations and labels
ax.set_xticks(xticks, crs=ccrs.PlateCarree()) # set ticks
ax.set_xticklabels(xticklabels) # set tick labels
# place y-ticks
yticks,yticklabels = ocean_row.tick_single_ax(ylims,ytickspacing) # tick locations and labels
ax.set_yticks(yticks, crs=ccrs.PlateCarree()) # set ticks
ax.set_yticklabels(yticklabels) # set tick labels
return ax
@staticmethod
def tick_single_ax(lims,tickstep,dropticksmid = False):
'''
Create ticks and ticklabels for tick_placer() for a single ax
Parameters
----------
lims: list [2]
axis limits
tickstep: scalar
tick spacing [deg]
dropticksmid: boolean
drop middle ticks when True
Returns
-------
ticks: vector
tick locations [deg]
ticklables: vector
ticklabels [deg]
'''
ticks = np.arange(np.ceil(lims[0]/tickstep)\
*tickstep,np.floor((lims[1])/tickstep+1)\
*tickstep,tickstep) # range of ticks
tickstepstr = str(tickstep) # ticks to dring
ticksdec = tickstepstr[::-1].find('.') # find decimal numbers
if ticksdec < 0: # if no decimal numbers
ticksdec = 0 # zero decimal numbers
ticklabels = [("%."+str(ticksdec)+"f"+"$^\circ$") % number for number in ticks] # make tick labels
if dropticksmid: # when mid ticks are dropped
ticklabels[1:-1] = [' '] * len(ticklabels[1:-1]) # drop mid ticks
return ticks,ticklabels # place ticks along axis
@staticmethod
def plt_base_simple(xlims,ylims,ytickspacing,xtickspacing):
'''
Create base of plot
Parameters
----------
xlims: list [2]
x-axis limits [deg]
xlims: list [2]
y-axis limits [deg]
xtickspacing: scalar
x tick spacing [deg]
ytickspacing: scalar
y tick spacing [deg]
Returns
-------
fig: fig
created figure
ax: ax
axes of created figure
'''
fig = plt.figure() # create figure
plateCr = ccrs.PlateCarree()
plateCr._threshold = plateCr._threshold/10. #set finer threshold
ax = plt.axes(projection=plateCr)
ax.set_extent([xlims[0],xlims[1],ylims[0],ylims[1]], ccrs.PlateCarree()) # set limits
ax.add_feature(cfeature.LAND,facecolor=[0.95,0.95,0.95]) # add land
ax.add_feature(cfeature.COASTLINE) # add coast
ax = ocean_row.tick_placer(ax,xlims,ylims,xtickspacing,ytickspacing) # place ticks
return fig,ax
@staticmethod
def plt_hindcast(hc,ptype):
'''
Plot montly average of hindcast data
Parameters
----------
hc: dictonary
contains all hindcast datasets
dsa: air
dsw: water
dswh: wave height
ptype: string
plot type
water velocity
water temperature
air velocity
daily wave height max
daily wave height mean
Returns
-------
fig: fig
created figure
ax: ax
axes of created figure
'''
tokts = 1.94384449 # convert m/s to kts
xlims = [30,130] # longitude limits [deg]
ylims = [-40,0] # latitude limits [deg]
if ptype == 'water velocity':
ds = hc['dsw']
months = pd.DatetimeIndex(ds.time.values).month
u = ds.uw*tokts
v = ds.vw*tokts
z = np.sqrt(u**2+v**2)
pmin = 0
pmax = 1.5
qscale = .2
unit = 'kts'
elif ptype == 'water temperature':
ds = hc['dsw']
months = pd.DatetimeIndex(ds.time.values).month
z = ds.Tw
pmin = 18
pmax = 32
unit = 'K'
elif ptype == 'air velocity':
ds = hc['dsa']
months = pd.DatetimeIndex(ds.time.values).month
u = ds.ua*tokts
v = ds.va*tokts
z = np.sqrt(u**2+v**2)
pmin = 0
pmax = 15
qscale = 2
unit = 'kts'
elif ptype == 'daily wave height max':
ds = hc['dswh']
months = pd.DatetimeIndex(ds.time.values).month
z = ds.H_day_max
pmin = 0
pmax = 4
unit = 'm'
elif ptype == 'daily wave height mean':
ds = hc['dswh']
months = pd.DatetimeIndex(ds.time.values).month
z = ds.H_day_mean
pmin = 0
pmax = 4
unit = 'm'
x = ds.lon.values
y = ds.lat.values
X,Y = np.meshgrid(x,y)
for i in range(1,13):
zi = z.isel(time=(months==i)).mean('time').values
fig,ax = ocean_row.plt_base_simple(xlims,ylims,10,10)
im = ocean_row.plt_modified_imshow(x,y,zi,ax=ax,cb=True,vmin=pmin,vmax=pmax,cmap='turbo')
if (ptype == 'water velocity') | (ptype == 'air velocity'):
ui = u.isel(time=(months==i)).mean('time').values
vi = v.isel(time=(months==i)).mean('time').values
qidx = int(np.ceil(X.shape[1]/30))
ax.quiver(X[::qidx,::qidx],Y[::qidx,::qidx],ui[::qidx,::qidx],vi[::qidx,::qidx],units='xy',scale_units='xy',scale=qscale)
ax.set(xlim=xlims,ylim=ylims,title=f'{calendar.month_name[i]}, {ptype} ({unit})')
plt.show()
pass
def plt_vrow_wind_only(self):
'''
Plot velocity depending on head and tail wind, for a given speed in perfect conditions (no wind and no waves)
Parameters
----------
self
Returns
-------
fig: fig
created figure
ax: ax
axes of created figure
'''
vperf = self.mcf['start_lon']
Gw = self.mcf['Gw']
GwA = self.mcf['GwA']
Ga = self.mcf['Ga']
Frmag = self.rowing_force() # calculate total rowing force
G = Ga/Gw # ratio of air friction cosntant to water friction constant
tokts = 1.94384449 # convert m/s to kts
Ngrid = 1000 # number of air velocities
vairr = np.linspace(-30,30,Ngrid)/tokts # grid of air velocities
# allocate
vrowairr = np.empty(Ngrid) # boat speed at vperf
for i in range(Ngrid): # loop over all grid points
# calculate boat speeds
vrowairr[i] = fsolve(ocean_row.total_force,\
list((vperf+vairr[i]*0.1)*np.array([1,0])),\
args=(np.array([0,0]),np.array([vairr[i],0]),np.array([1,0]),Frmag,True,False,Gw,GwA,Ga))[0]
fig,ax = plt.subplots(figsize=(4,6)) # create figure
ax.axhline(y=0, color='black', linestyle='-',linewidth=.5) # add line at y=0
ax.plot(np.abs(vairr*tokts),vrowairr*tokts,color='blue',linewidth=2) # plot rowing speed as function of wind speed
ax.set(xlim=[0,np.max(np.abs(vairr))*tokts],ylim=[-3,5],title='$G = {:.4f}$'.format(G),xlabel='Wind speed (kts)',ylabel='Rowing speed (kts)') # add labels
return fig,ax
@staticmethod
def angle_between(p1,p2):
'''
Calculate angle between two points
Parameters
----------
p1: vector [2]
point 1 [x,y]
p2: vector [2]
point 2 [x,y]
Returns
-------
angle: scalar
anlge between points
'''
ang1 = np.arctan2(*p1[::-1]) # take arctan of first point
ang2 = np.arctan2(*p2[::-1]) # take arctan of second point
angle = np.rad2deg((ang1 - ang2) % (2 * np.pi)) # calculate angle between points
return angle
@staticmethod
def filter_circular(Xa,Xt,Np,keepsides=True,geo=False):
'''
Filter points close to target (bin by anlge with target)
Parameters
----------
Xa: list
boat coordinates [2]
Xt: vector [2]
target coordinate
Np: scalar
number of remaining locations after filtering
keepsides: boolean
most outer locations are kept in addition when True
geo: boolean
True: when applied to coordinates (in degree)
False: when applied to locations (in m)
Returns
-------
Xaf: list
filtered locations [2]
idxkeep: boolean vector
indexes of Xa that remain (not removed by the filtration)
'''
NXa = Xa.shape[0] # number of points
if NXa>Np: # apply this function if there are more points then requested segments
Xa0 = Xa-Xt # set Xt as centre point (0,0)
if geo: # when applied to coordinates
dis = np.array([Geodesic.WGS84.Inverse(x[1],x[0],Xt[1],Xt[0])['s12'] for x in Xa]) # calculate distance to target
phi = np.array([Geodesic.WGS84.Inverse(Xt[1],Xt[0],x[1],x[0])['azi1']/180*np.pi for x in Xa]) # calculate angle from target to boat
else: # when applied to locations
dis = np.sqrt(Xa0[:,0]**2+Xa0[:,1]**2) # calculate distance to target
phi = np.arctan2(Xa0[:,1],Xa0[:,0])+np.pi # calculate angle from target to boat
# find outer locations by finding the interval using the largest difference between angles
phis = np.sort(phi) # sort angles
phis2 = np.append(phis,phis+2*np.pi) # add angles 2pi displaced
idxmdif = np.argmax(np.diff(phis2)[:NXa]) # find idx of largest angle
# define interval
if idxmdif == NXa-1: # when idx of largest angle is the last sorted location
# interval boundaries
phimin = phis[0] # start of interval
phimax = phis[idxmdif] # end of interval
# normalize interval by setting phimin to zero
phinmin = 0 # normalized start of interval
phinmax = phimax-phimin # normalized end of interval
else: # when idx of largest angle is in between the sorted locations
# interval boundaries
phimin = phis[idxmdif+1] # start of interval
phimax = phis[idxmdif] # end of interval
# normalize interval by setting phimin to zero
phinmin = 0 # normalized start of interval
phinmax = 2*np.pi-(phimin-phimax) # normalized end of interval
# normalize all angles
phin = phi-phimin # normalize all angles
phin[phin<0] += 2*np.pi # updated to boundaries
phin[phin>2*np.pi] -= 2*np.pi # updated to boundaries
Nseg = Np # number of segments
Nsegd = 0 # allocate number of filled segments (should equal Np)
flag1 = True
flag2 = True
while (flag1 | flag2): # when there are less filled segments than requested number of locations
binedge = np.linspace(0,phinmax,Nseg+1) # bin segments
inds = np.digitize(phin,binedge)-1 # check if locoations in segments
inds[inds>(Nseg-1)] = Nseg-1 # include angle at phimax in last segment
indsunq = np.unique(inds) # find unique segments
Nsegd = len(indsunq) # total number of filled segments
# if there are less filled segments (which meanch less remainging locations) than requested locations
if flag1:
if (Nsegd < Np):
Nseg += 1 # increase number of segments
else:
flag1 = False
if flag1 == False:
# allocate array to indicate if location should be kept (= True) or removed (= False)
idxkeep = np.zeros(NXa).astype(bool)
idxside = np.zeros(NXa).astype(bool)
idxclose = np.zeros(NXa).astype(bool)
for i in range(Nsegd): # loop over all segments
idx = inds == indsunq[i] # get indexes of segment
Xa0i = Xa0[idx,:] # coordinates of locations in segment
disi = dis[idx] # distances of locations in segment
idxmin = np.argmin(disi) # find which location is closest to target
idxclose[np.where(idx)[0][idxmin]] = True # keep closest location
if keepsides: # when keepsides is requested
idxside[np.argmax(phin)] = True # keep minumum anlge
idxside[np.argmin(phin)] = True # keep maximum angle
idxkeep = idxclose | idxside
else:
idxkeep = idxclose
if np.sum(idxkeep)>Np:
Nseg -= 1
else:
flag2 = False
Xa0b = Xa0[idxkeep] # filter locations
Xaf = Xa0b+Xt # transfer locations (with respect to target) back to original locations
else: # if there are less available points then requested segments, then no filtering is required
Xaf = Xa # keep all locations
idxkeep = np.ones(NXa).astype(bool) # note all indexes as kept
return Xaf,idxkeep
@staticmethod
def get_hindcast_at_loc(X_a,datet,dsa,dsw,dsa_time,dsa_lon,dsa_lat,dsw_time,dsw_lon,dsw_lat):
'''
Get the hindcast velocities (air and water) at given boat locations
Parameters
----------
X_a: matrix
N boat locations [N,2]
datet: datetime
time
dsa: xarray
full dataset air
dsw: xarray
full dataset water
dsa_time: vector
coordinate
dsa_lon: vector
coordinate
dsa_lat: vector
coordinate
dsw_time: vector
coordinate
dsw_lon: vector
coordinate
dsw_lat: vector
coordinate
Returns
-------
Va_a: matrix
air velocity at N boat locations [N,2]
Vw_a: matrix
water velocity at N boat locations [N,2]
'''
N = X_a.shape[0]
iatime = np.argmin(np.abs(dsa_time-datet))
iwtime = np.argmin(np.abs(dsw_time-datet))
Va_a = np.empty([N,2])
Vw_a = np.empty([N,2])
for iX in range(N): # loop over all locations
Xi = X_a[iX,:]
ialon = np.argmin(np.abs(dsa_lon-Xi[0]))
ialat = np.argmin(np.abs(dsa_lat-Xi[1]))
iwlon = np.argmin(np.abs(dsw_lon-Xi[0]))
iwlat = np.argmin(np.abs(dsw_lat-Xi[1]))
Va_a[iX,0] = dsa.isel(time=iatime,lon=ialon,lat=ialat).ua.item()
Va_a[iX,1] = dsa.isel(time=iatime,lon=ialon,lat=ialat).va.item()
Vw_a[iX,0] = dsw.isel(time=iwtime,lon=iwlon,lat=iwlat).uw.item()
Vw_a[iX,1] = dsw.isel(time=iwtime,lon=iwlon,lat=iwlat).vw.item()
return Va_a,Vw_a
def run_ocean_row(self,hc,dmode,imax,rtmin):
'''
Run the ocean rowing model
Parameters
----------
self
hc: dictonary
contains all hindcast datasets
dsa: air
dsw: water
dswh: wave height
dmode: dictionary
mode: string
rowing mode (destination, track, optimized)
if mode == destination
anchor_drop: boolean
True: always drop anchor when resting
False: never drop anchor when resting
filters_on: boolean
True: remove location when on land
False: no filtering based on land
if mode == track
Xa_pr: matrix [N_pr,2]
Locations of track
tod_0: string *OPTIONAL*
string of datetime to start row
overwrites tod_0 of mcf in function
anchor_drop: boolean
True: always drop anchor when resting
False: never drop anchor when resting
filters_on: boolean
True: remove location when on land
False: no filtering based on land
if mode == optimze
angle_max: scalar
maximum rowing angle w.r.t. destination
angle_Nstep: scalar
number of rowing directions (should be odd)
Np: scalar
number of filtered points
imax: scalar
maximum number of iterations
rtmin: scalar
threshold of distance till target [m]
Returns
-------
Xs: matrix [Nend,Np,2]
location [deg]
Vas: matrix [Nend,Np,2]
velocity air
Vws: matrix [Nend,Np,2]
velocity water
Vs: matrix [Nend,Np,2]
velocity boat (during previous segment)
datets: vector [Nend]
datetime
Nhists: matrix [Nend,Np]
history of previous location
Ns: vector [Nend]
number of boat locations per iteration
'''
# values from self
vperf = self.mcf['vperf']
Gw = self.mcf['Gw']
GwA = self.mcf['GwA']
Ga = self.mcf['Ga']
tod_start = self.mcf['tod_start']
tod_stop = self.mcf['tod_stop']
pdom = ocean_row.create_domain() # create allowed domain to row
mode = dmode['mode'] # get mode
if mode == 'destination':
Np = 1 # max number of points after filtering
anchor_drop = dmode['anchor_drop'] # configure anchor
filters_on = dmode['filters_on']
elif mode == 'track':
Np = 1 # max number of points after filtering
anchor_drop = dmode['anchor_drop'] # configure anchor
Xa_pr = dmode['Xa_pr'] # all previous row locations [deg]
N_pr = Xa_pr.shape[0] # total number of locations
filters_on = dmode['filters_on']
elif mode == 'optimize':
tod_0 = self.mcf['tod_0'] # starting datetime
phir = np.linspace(-dmode['angle_max'],dmode['angle_max'],dmode['angle_Nstep']) # angle offsets to target
Nphir = len(phir) # number of angle offsets
Np = dmode['Np'] # max number of points after filtering
filters_on = True
if mode == 'track':
if 'tod_0' in dmode: # if starting datetime is listed in dmode
tod_0 = dmode['tod_0'] # starting datetime
else: # take starting datetime from mcf
tod_0 = self.mcf['tod_0'] # starting datetime
start_lon = Xa_pr[0,0] # starting longitude
start_lat = Xa_pr[0,1] # starting latitude
stop_lon = Xa_pr[N_pr-1,0] # destination longitude
stop_lat = Xa_pr[N_pr-1,1] # destination latitude
else:
tod_0 = self.mcf['tod_0'] # starting datetime
start_lon = self.mcf['start_lon'] # starting longitude
start_lat = self.mcf['start_lat'] # starting latitude
stop_lon = self.mcf['stop_lon'] # destination longitude
stop_lat = self.mcf['stop_lat'] # destination latitude
# intialize model
# calculate force from velocity in perfect conditions
Frmag = self.rowing_force() # calculate total rowing force
int0,introw,intdur,intN,k,dt = self.rowing_activity() # row activity
datet0 = np.array(np.datetime64(tod_0,'s'))
# distances
X0 = np.array([start_lon,start_lat]) # starting point
Xt = np.array([stop_lon,stop_lat]) # target point
rt0 = Geodesic.WGS84.Inverse(X0[1],X0[0],Xt[1],Xt[0])['s12'] # distance to target
# intialize
N = 1 # number of starting points
No = 1 # number of options
X_a = np.empty([N,2]) # allocate all locations
X_a[0,:] = X0 # location
i = 0 # iteration
rto = np.array([Geodesic.WGS84.Inverse(x[1],x[0],Xt[1],Xt[0])['s12'] for x in X_a]) # geodesic distance to target
datet = datet0 # datetime
Nhist = np.nan # previous boat location
V_a = np.array([0,0]) # boat velocity
# datasets
dsa = hc['dsa'] # hindcast air speeds
dsw = hc['dsw'] # hindcast water speed
# dataset coordinates
dsa_time = dsa.time.values
dsa_lon = dsa.lon.values
dsa_lat = dsa.lat.values
dsw_time = dsw.time.values
dsw_lon = dsw.lon.values
dsw_lat = dsw.lat.values
# allocate
Xs = np.empty([imax,Np,2]); Xs.fill(np.nan) # locations
Vas = np.empty([imax,Np,2]); Vas.fill(np.nan) # air velocity
Vws = np.empty([imax,Np,2]); Vws.fill(np.nan) # water velocity
Vs = np.empty([imax,Np,2]); Vs.fill(np.nan) # boat velocity
datets = np.empty([imax],dtype='datetime64[s]'); # time series
Nhists = np.zeros([imax,Np]) # history of previous boat location
Ns = np.zeros([imax]) # number of boat locations
# create plot
xlims = [30,130] # longitude limits [deg]
ylims = [-40,0] # latitude limits [deg]
fig,ax = ocean_row.plt_base_simple(xlims,ylims,10,10) # create figure
plt.gcf().set_size_inches(13,6)
ax.set_title(pd.to_timedelta(np.round((datet-datet0)/np.timedelta64(1, 's')),'s'))
circt = np.array([[Geodesic.WGS84.Direct(Xt[1],Xt[0],angle,rtmin)[key] for key in ['lon2','lat2']] for angle in np.arange(0,360,1)])
ax.plot(circt[:,0],circt[:,1],'r-')
ax.plot(X0[0],X0[1],'o',color='blue',markersize=7)
ax.plot(Xt[0],Xt[1],'o',color='red' ,markersize=7)
if mode == 'optimize':
pltraj = ax.plot(X_a[:,0],X_a[:,1],'o',color='r',markersize=2)
else:
pltraj = ax.plot(Xs[:i+1,0,0],Xs[:i+1,0,1],'r-',linewidth=1)
if dmode['mode']=='track':
plt.plot(Xa_pr[:,0],Xa_pr[:,1],'k-',linewidth=1)
while (i < (imax-1)) & np.all(rto > rtmin): # start loop
# get hindcast data
Va_a,Vw_a = self.get_hindcast_at_loc(X_a,datet,dsa,dsw,dsa_time,dsa_lon,dsa_lat,dsw_time,dsw_lon,dsw_lat)
# set hindcast velocities to zero when no data is available
Va_a[np.isnan(Va_a)] = 0
Vw_a[np.isnan(Vw_a)] = 0
# store values
Xs [i,0:N,:] = X_a
Vas [i,0:N,:] = Va_a
Vws [i,0:N,:] = Vw_a
Vs [i,0:N,:] = V_a
Nhists[i,0:N] = Nhist
Ns [i] = N
datets[i] = datet
i += 1 # next iteration
# rowing activity
rowt = introw[k] # check rowing activity options
if (mode == 'destination')|(mode == 'track'):
No = 1 # number of options
if rowt == True: # if rowing
phio = np.array([0]) # rowing direction options
rowo = np.array([True]) # rowing activity
anco = np.array([False]) # achor
if rowt == False: # if resting
phio = np.array([np.nan]) # rowing direction options
rowo = np.array([False]) # rowing activity
anco = np.array([anchor_drop]) # anchor
elif mode == 'optimize':
if rowt == True: # if rowing
No = Nphir+2 # can row (Nphir options) or rest(2 options: with and without anchor)