Particles simulating the oil spill in English Bay, Vancouver on April 8, 2015 or April 9 in UTC.
In [1]:
from IPython.display import Image
import netCDF4 as nc
import numpy as np
from salishsea_tools import (nc_tools, viz_tools, tidetools)
import matplotlib.pyplot as plt
import matplotlib.cm as cm
%matplotlib inline
In [2]:
Image(filename='/ocean/imachuca/MEOPAR/Ariane/results/English_Bay_Spill/Oil_Spill1.jpg')
Out[2]:
In [3]:
Image(filename='/ocean/imachuca/MEOPAR/Ariane/results/English_Bay_Spill/Oil_Spill2.png')
Out[3]:
In [4]:
grid = nc.Dataset('../../NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
day1U = nc.Dataset('/data/dlatorne/MEOPAR/SalishSea/nowcast/09apr15/SalishSea_1h_20150409_20150409_grid_U.nc','r')
day2U = nc.Dataset('/data/dlatorne/MEOPAR/SalishSea/nowcast/10apr15/SalishSea_1h_20150410_20150410_grid_U.nc','r')
day4U = nc.Dataset('/data/dlatorne/MEOPAR/SalishSea/nowcast/12apr15/SalishSea_1h_20150412_20150412_grid_U.nc','r')
day1V = nc.Dataset('/data/dlatorne/MEOPAR/SalishSea/nowcast/09apr15/SalishSea_1h_20150409_20150409_grid_V.nc','r')
ugrid = day1U.variables['vozocrtx']
vgrid = day1V.variables['vomecrty']
zlevels = day1U.variables['depthu']
timesteps = day1U.variables['time_counter']
t=0
timestamp1 = nc_tools.timestamp(day1U, t)
timestamp2 = nc_tools.timestamp(day2U, t)
timestamp4 = nc_tools.timestamp(day4U, t)
In [5]:
bathy, X, Y = tidetools.get_bathy_data(grid)
lat = 49.289227
lon = -123.187080
yind, xind = tidetools.find_closest_model_point(lon, lat, X, Y, bathy)
print xind, yind
333 455
In [6]:
zlevel = 0
x_slice = np.arange(260, 370)
y_slice = np.arange(420, 520)
ugrid_tzyx = np.ma.masked_values(ugrid[t, zlevel, y_slice, x_slice], 0)
vgrid_tzyx = np.ma.masked_values(vgrid[t, zlevel, y_slice, x_slice], 0)
u_tzyx, v_tzyx = viz_tools.unstagger(ugrid_tzyx, vgrid_tzyx)
fig, ax = plt.subplots(1, 1, figsize=(20, 10))
viz_tools.set_aspect(ax)
ax.streamplot(x_slice[1:], y_slice[1:], u_tzyx, v_tzyx, density=2)
ax.plot(xind, yind, marker='o', markersize=14, markeredgewidth=2, color='red', label='Origin')
viz_tools.plot_land_mask(ax, grid, xslice=x_slice, yslice=y_slice)
ax.set_xlim(x_slice[0], x_slice[-1])
ax.set_ylim(y_slice[0], y_slice[-1])
ax.set_xlabel('x Index')
ax.set_ylabel('y Index')
timestamp = nc_tools.timestamp(day1U, t)
ax.set_title('Velocity Field\n'+timestamp1.format('DD-MMM-YYYY, HH:mm [UTC]'))
legend = ax.legend(numpoints=1, loc=4)
Ariane: delta_t=60, freq=15, tunit=3600, ntfic=1, lmt=24, nb_output=48, fl=0.5, 4.5, 8.5, nmax=3
In [7]:
part1 = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/English_Bay_Spill/0_4_8_12hrs/ariane_trajectories_qualitative.nc','r');
lon1=part1.variables['traj_lon']
lat1=part1.variables['traj_lat']
dep1=part1.variables['traj_depth']
tim1=part1.variables['traj_time']
inx1=part1.variables['init_x']
iny1=part1.variables['init_y']
int1=part1.variables['init_t']
part2 = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/English_Bay_Spill/12_16_20_12hrs/ariane_trajectories_qualitative.nc','r');
lon2=part2.variables['traj_lon']
lat2=part2.variables['traj_lat']
dep2=part2.variables['traj_depth']
tim2=part2.variables['traj_time']
inx2=part2.variables['init_x']
iny2=part2.variables['init_y']
int2=part2.variables['init_t']
part3 = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/English_Bay_Spill/0_48hrs/ariane_trajectories_qualitative.nc','r');
lon3=part3.variables['traj_lon']
lat3=part3.variables['traj_lat']
dep3=part3.variables['traj_depth']
tim3=part3.variables['traj_time']
inx3=part3.variables['init_x']
iny3=part3.variables['init_y']
int3=part3.variables['init_t']
partt = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/English_Bay_Spill/test_15min/ariane_trajectories_qualitative.nc','r');
lont=partt.variables['traj_lon']
latt=partt.variables['traj_lat']
dept=partt.variables['traj_depth']
timt=partt.variables['traj_time']
inxt=partt.variables['init_x']
inyt=partt.variables['init_y']
intt=partt.variables['init_t']
part4 = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/English_Bay_Spill/48hrs_multiples/ariane_trajectories_qualitative.nc','r');
lon4=part4.variables['traj_lon']
lat4=part4.variables['traj_lat']
dep4=part4.variables['traj_depth']
tim4=part4.variables['traj_time']
inx4=part4.variables['init_x']
iny4=part4.variables['init_y']
int4=part4.variables['init_t']
In [8]:
def plot_map(ax, grid):
fig, ax = plt.subplots(1,1,figsize=(20,15))
viz_tools.plot_land_mask(ax,grid,coords='map')
viz_tools.set_aspect(ax)
ax.set_xlim([-123.29,-123.13])
ax.set_ylim([49.22,49.34])
return ax
In [32]:
def add_labels(ax):
ax.plot(-123.155, 49.274, marker='o', markersize=14, markeredgewidth=2, color='MistyRose') #Kitsilano Beach
ax.plot(-123.15, 49.30, marker='o', markersize=14, markeredgewidth=2, color='MistyRose') #Third Beach
ax.plot(-123.216, 49.277, marker='o', markersize=14, markeredgewidth=2, color='MistyRose') #Spanish Banks
bbox_args = dict(boxstyle='square', facecolor='white')
ax.annotate('Spanish\nBanks', (-123.22, 49.265), fontsize=15, color='white')
ax.annotate('Kitsilano\nBeach', (-123.155, 49.26), fontsize=15, color='white')
ax.annotate('Third\nBeach', (-123.145, 49.292), fontsize=15, color='white')
return ax
In [35]:
ax = plot_map(ax, grid)
ax = add_labels(ax)
ax.set_title('Particles originating at point of oil spill in English Bay\n'+timestamp1.format('DD-MMM-YYYY')
+', Duration = 12 hrs, Frequency = 15 min' )
legend.get_title().set_fontsize('20')
ax.annotate('Origin', (-123.19, 49.282), fontsize=20)
n = np.arange(lon1.shape[1]+lon2.shape[1])
colors = cm.rainbow(np.linspace(0, 1, len(n)))
for N, c in zip(n[0:lon1.shape[1]], colors[0:3]):
ax.plot(lon1[1:,N],lat1[1:,N], linewidth=3.0, color=c, label=int1[N]-0.5)
ax.plot(lon1[0,N],lat1[0,N], marker='o', markersize=14, markeredgewidth=2, color='MistyRose')
for N, c in zip(n[0:lon2.shape[1]], colors[3:7]):
ax.plot(lon2[1:,N],lat2[1:,N], linewidth=3.0, color=c, label=int2[N]-0.5)
ax.plot(lon2[0,N],lat2[0,N], marker='o', markersize=14, markeredgewidth=2, color='MistyRose')
ax.legend(loc=4, ncol=3, title=r'Initial Times of Particle Trajectories [UTC]')
Out[35]:
<matplotlib.legend.Legend at 0x7f7cd9dace50>
In [34]:
ax = plot_map(ax, grid)
ax = add_labels(ax)
ax.set_title('Particles originating at point of oil spill in English Bay\n'+timestamp1.format('DD-MMM-YYYY')+' to '+timestamp2.format('DD-MMM-YYYY')
+', Duration = 48 hrs, Frequency = 15 min' )
ax.plot(lon3[1:,0],lat3[1:,0], linewidth=3.0, color='red')
ax.plot(lon3[0,0],lat3[0,0], marker='o', markersize=14, markeredgewidth=2, color='MistyRose', label='Origin')
ax.annotate('Origin', (-123.19, 49.282), fontsize=20)
Out[34]:
<matplotlib.text.Annotation at 0x7f7cd9db8b10>
In [12]:
times = np.arange(0.5,48.75,0.25)
len_times = len(times)
times_ind = np.arange(len_times)
times_ind_r = times_ind[::-1]
times_ind_r
Out[12]:
array([192, 191, 190, 189, 188, 187, 186, 185, 184, 183, 182, 181, 180,
179, 178, 177, 176, 175, 174, 173, 172, 171, 170, 169, 168, 167,
166, 165, 164, 163, 162, 161, 160, 159, 158, 157, 156, 155, 154,
153, 152, 151, 150, 149, 148, 147, 146, 145, 144, 143, 142, 141,
140, 139, 138, 137, 136, 135, 134, 133, 132, 131, 130, 129, 128,
127, 126, 125, 124, 123, 122, 121, 120, 119, 118, 117, 116, 115,
114, 113, 112, 111, 110, 109, 108, 107, 106, 105, 104, 103, 102,
101, 100, 99, 98, 97, 96, 95, 94, 93, 92, 91, 90, 89,
88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76,
75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63,
62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50,
49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37,
36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24,
23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11,
10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0])
Ariane: delta_t=60, freq=15, tunit=3600, ntfic=1, lmt=120, nb_output=192 (april 9-13), fl=0.5 to 48.5, nmax=193
file = open("initial_positions_write.txt", "w") for n, time in zip(times_ind, times): file.write("333 455 -1.0 " +str(times[n])+ " 1.0\n") file.close()
In [36]:
ax = plot_map(ax, grid)
ax = add_labels(ax)
ax.annotate('Origin', (-123.19, 49.282), fontsize=20)
ax.set_title('Particles originating at point of oil spill in English Bay\n'+
'One particle dropped every 15 minutes from '
+timestamp1.format('DD-MMM-YYYY') + ', 00:00'
+' to '+timestamp2.format('DD-MMM-YYYY') + ', 24:00\n'
+'Duration = Trajectories end on '+timestamp2.format('DD-MMM-YYYY') + ', 24:00'
+', Frequency of particle positions = 15 min' )
for i, j in zip(times_ind, times_ind_r):
ax.plot(lon4[0:j+1,i],lat4[0:j+1,i], color='DodgerBlue')
ax.plot(lon4[0,0], lat4[0,0], marker='o', markersize=14, markeredgewidth=2, color='MistyRose')
plt.show()
In [37]:
ax = plot_map(ax, grid)
ax = add_labels(ax)
ax.annotate('Origin', (-123.19, 49.282), fontsize=20)
ax.set_title('Particles originating at point of oil spill in English Bay\n'+
'One particle dropped every 15 minutes from '
+timestamp1.format('DD-MMM-YYYY')+', 00:00'
+' to '+timestamp2.format('DD-MMM-YYYY')+', 24:00\n'
+'Duration = 48 hrs (Final trajectory ends on '+timestamp4.format('DD-MMM-YYYY')+', 24:00\n'
+'Frequency of particle positions = 15 min' )
ax.plot(lon4[1:,:],lat4[1:,:], color='DodgerBlue')
ax.plot(lon4[0,0], lat4[0,0], marker='o', markersize=14, markeredgewidth=2, color='MistyRose')
plt.show()
In [63]:
ax = plot_map(ax, grid)
ax = add_labels(ax)
ax.annotate('Origin', (-123.19, 49.282), fontsize=20)
legend.get_title().set_fontsize('20')
ax.set_title('Particles originating at point of oil spill in English Bay\n'+
'One particle dropped every 15 minutes from '
+timestamp1.format('DD-MMM-YYYY')+', 00:00'
+' to '+timestamp2.format('DD-MMM-YYYY')+', 24:00\n'
+'Duration = 48 hrs (Final trajectory ends on '+timestamp4.format('DD-MMM-YYYY')+', 24:00\n'
+'Frequency of particle positions = 15 min' )
ax.plot(lon4[1:, 0],lat4[1:, 0], color='MediumSlateBlue', label='09-Apr')
ax.plot(lon4[1:, 0:96],lat4[1:, 0:96], color='MediumSlateBlue')
ax.plot(lon4[1:, 96],lat4[1:, 96], color='DodgerBlue', label='10-Apr')
ax.plot(lon4[1:, 96:],lat4[1:, 96:], color='DodgerBlue')
ax.plot(lon4[0,0], lat4[0,0], marker='o', markersize=14, markeredgewidth=2, color='MistyRose')
ax.legend(loc=4, ncol=1, title=r'Release date')
plt.show()
In [ ]: