-
Notifications
You must be signed in to change notification settings - Fork 0
/
1_plot_scenario_maps.py
198 lines (175 loc) · 8.94 KB
/
1_plot_scenario_maps.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
'''
Script to plot mapos of release locations for all 3 Parcels scenarios in one figure
Author: vesnaber
Additinal files needed: 1) croco_grd.nc - input grid file
from CROCO model (see repository of SCARIBOS)
2) Release locations, generated by the following scripts:
- 1_release_locations_COASTCON.py
- 1_release_locations_REGIOCON.py
- 1_release_locations_HOTSPOTS.py
3) Buffered areas for Intra-island connectivity scenario (i.e. COASTCON),
named COASTCON_buffered_areas.npy, generated by:
- 1_buffer_areas_COASTCON.py
'''
# Import libraries
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Point, MultiPoint
import matplotlib.ticker as ticker
from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
# Upload grid from model CROCO (see repository of SCARIBOS) - for bathymetry
config = 'SCARIBOS_V8' # CROCO model configuration name
path = '~/croco/CONFIG/' + config + '/CROCO_FILES/'
grid = xr.open_dataset(path + 'croco_grd.nc')
bathymetry = grid.h.values
landmask = grid.mask_rho.values
landmask = np.where(landmask == 0, 1, np.nan)
# COASTCON (Intra-island connectivity):
buffered_areas = np.load('COASTCON/buffer_areas/COASTCON_buffered_areas.npy', allow_pickle=True)
release_locations = ['zone_1', 'zone_2', 'zone_3', 'zone_4', 'zone_5', 'zone_6', 'zone_7', 'zone_8']
colors = plt.cm.ocean_r(np.linspace(0, 1, len(release_locations)+1))
# HOTSPOTS:
hotspots = np.load('INPUT/release_locs_HOTSPOTS.npy')
xmin, xmax, ymin, ymax = -69.49, -68.49, 11.67, 12.67 # limits for plotting hotspots
# REGIOCON (Coastal connectivity):
ARUB = np.load('INPUT/release_locs_REGIOCON_ARUB.npy')
BONA = np.load('INPUT/release_locs_REGIOCON_BONA.npy')
VEIS = np.load('INPUT/release_locs_REGIOCON_VEIS.npy')
VECO = np.load('INPUT/release_locs_REGIOCON_VECO.npy')
curacao_coast_points = np.load('INPUT/coastal_points_CURACAO.npy')
# make buffer area from coastal points of Curaçao
points = [Point(lon, lat) for lon, lat in zip(curacao_coast_points[0], curacao_coast_points[1])]
multi_point = MultiPoint(points)
buffer_distance = 0.05 # Define the buffer distance
buffered_area_REGIO = multi_point.buffer(buffer_distance)
buffered_gdf = gpd.GeoDataFrame(geometry=[buffered_area_REGIO])
# Create a polygon from the shapely geometry of buffered_area_REGIO
polygon = Polygon(list(zip(buffered_area_REGIO.exterior.xy[0], buffered_area_REGIO.exterior.xy[1])),
facecolor='yellow', edgecolor='orange', hatch='///', alpha=0.6,linewidth=1, label='Destination area: Curaçao', zorder=1)
# figure settings:
font0 = 14
font1 = 13
font2 = 16
cura_color = 'slateblue'
# bathymetry colorscheme
def custom_div_cmap(numcolors=50, name='custom_div_cmap',
mincol='blue', midcol2='yellow', midcol='white', maxcol='red'):
""" Create a custom diverging colormap with three colors
Default is blue to white to red with 11 colors. Colors can be specified
in any way understandable by matplotlib.colors.ColorConverter.to_rgb()
"""
cmap = LinearSegmentedColormap.from_list(name=name,
colors=[mincol, midcol, midcol2, maxcol],
N=numcolors)
return cmap
blevels = [-5000, -4000, -3000, -2500, -2000, -1500, -1000,-800, -600, -400, -200, 0] # define levels for plotting (transition of colorbar)
N = (len(blevels)-1)*2
bathy_cmap = custom_div_cmap(N, mincol='#696969', midcol='dimgrey', midcol2='#888888' ,maxcol='w')
vmin = -8000
vmax = 0
# FIGURE
fig = plt.figure(figsize=(14, 14))
gs = fig.add_gridspec(3, 4, height_ratios=[1, 1.3, 0.1], width_ratios=[0.01, 1, 1,0.01])
# Top right plot: HOTSPOTS
ax2 = fig.add_subplot(gs[0, 1])
bathy = ax2.contourf(grid.lon_rho, grid.lat_rho, -bathymetry, 60, cmap=bathy_cmap, vmin=-8000, vmax=0)
for c in bathy.collections:
c.set_rasterized(True)
ax2.scatter(hotspots[0], hotspots[1], 0.9, color='darkblue', label='Release locations')
ax2.contourf(grid.lon_rho, grid.lat_rho, landmask, cmap='Greys', alpha=1, zorder=1)
ax2.contourf(grid.lon_rho, grid.lat_rho, landmask, cmap='Greys_r', alpha=0.9, zorder=2)
ax2.set_xlim(xmin-0.3, xmax+0.3)
ax2.set_ylim(ymin-0.3, ymax+0.3)
ax2.set_title('(a) Scenario 1: Hotspots', fontsize=font2)
ax2.set_aspect('equal', 'box')
ax2.set_xticks(np.arange(-69.5, -68, 0.5))
ax2.set_yticks(np.arange(11.5, 13.5, 0.5))
ax2.tick_params(axis='both', which='major', labelsize=font0)
ax2.legend(loc='upper right', fontsize=font1, markerscale=10)
ax2.set_xticklabels(['{:.1f}° W'.format(abs(x)) for x in ax2.get_xticks()])
ax2.set_yticklabels(['{:.1f}° N'.format(abs(y)) for y in ax2.get_yticks()])
# Top left plot: COASTCON release locations and destination zones
ax1 = fig.add_subplot(gs[0,2])
bathy = ax1.contourf(grid.lon_rho, grid.lat_rho, -bathymetry, 60, cmap=bathy_cmap, vmin=-8000, vmax=0)
for c in bathy.collections:
c.set_rasterized(True)
for idx, buffered_area in enumerate(buffered_areas):
ax1.fill(*buffered_area.exterior.xy, color=colors[idx+1], alpha=0.8)
ax1.contourf(grid.lon_rho, grid.lat_rho, landmask, cmap='Greys', alpha=1, zorder=1)
ax1.contourf(grid.lon_rho, grid.lat_rho, landmask, cmap='Greys_r', alpha=0.9, zorder=2)
ax1.set_xlim(xmin+0.15+0.1, xmax-0.15+0.1)
ax1.set_ylim(ymin+0.15+0.1, ymax-0.15+0.1)
ax1.set_title('(b) Scenario 2: Intra-island connectivity', fontsize=font2)
no_of_particles = []
for idx, release_location in enumerate(release_locations):
zone_data = np.load(f'INPUT/release_locs_COASTCON_{release_location}.npy')
ax1.scatter(zone_data[0], zone_data[1], 1, color='k')
no_of_particles.append(f'{len(zone_data[0])}')
ax1.legend([f'1: Klein Curaçao ({no_of_particles[0]} p.)', f'2: Oostpunt ({no_of_particles[1]} p.)',
f'3: Caracasbaai ({no_of_particles[2]} p.)', f'4: Willemstad ({no_of_particles[3]} p.)',
f'5: Bullenbaai ({no_of_particles[4]} p.)', f'6: Valentijnsbaai ({no_of_particles[5]} p.)',
f'7: Westpunt ({no_of_particles[6]} p.)', f'8: North Shore ({no_of_particles[7]} p.)'],
loc='upper right', fontsize=font1, markerscale=10)
ax1.set_aspect('equal', 'box')
ax1.set_xticks(np.arange(-69.2, -68.6, 0.2))
ax1.set_yticks(np.arange(12, 12.7, 0.2))
ax1.tick_params(axis='both', which='major', labelsize=font0)
ax1.set_xticklabels(['{:.1f}° W'.format(abs(x)) for x in ax1.get_xticks()])
ax1.set_yticklabels(['{:.1f}° N'.format(abs(y)) for y in ax1.get_yticks()])
# Bottom plot: REGIOCON release locations
ax3 = fig.add_subplot(gs[1, 1:])
bathy = ax3.contourf(grid.lon_rho, grid.lat_rho, -bathymetry, 60, cmap=bathy_cmap, vmin=-8000, vmax=0, extend='min')
for c in bathy.collections:
c.set_rasterized(True)
ax3.add_patch(polygon)
dotscoastsize = 4
ax3.contourf(grid.lon_rho, grid.lat_rho, landmask, cmap='Greys', alpha=1, zorder=1)
ax3.contourf(grid.lon_rho, grid.lat_rho, landmask, cmap='Greys_r', alpha=0.9, zorder=2)
ax3.scatter(ARUB[0], ARUB[1], dotscoastsize, color='orange', label='Aruba')
ax3.scatter(BONA[0], BONA[1], dotscoastsize, color='tomato', label='Bonaire')
ax3.scatter(VEIS[0], VEIS[1], dotscoastsize, color='dodgerblue', label='Venezuelan Islands')
ax3.scatter(VECO[0], VECO[1], dotscoastsize, color='lightseagreen', label='Venezuela (mainland)')
ax3.set_title('(c) Scenario 3: Coastal connectivity', fontsize=font2)
ax3.set_ylim(10.25, 13.5)
ax3.set_xlim(-70.5, -66)
ax3.legend(loc='lower left', fontsize=font1, markerscale=6)
ax3.set_xticks(np.arange(-70, -65, 1))
ax3.set_yticks(np.arange(11, 14, 1))
ax3.tick_params(axis='both', which='major', labelsize=font0)
ax3.xaxis.set_major_formatter(ticker.FormatStrFormatter('%.1f'))
ax3.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.1f'))
ax3.set_xticklabels(['{:.1f}° W'.format(abs(x)) for x in ax3.get_xticks()])
ax3.set_yticklabels(['{:.1f}° N'.format(abs(y)) for y in ax3.get_yticks()])
# colorbar
cbar = plt.colorbar(bathy, ax=ax3, orientation='vertical', pad=0.025, aspect=18, extend='min')
cbar.set_label('Depth [m]', fontsize=font1)
cbar.ax.tick_params(labelsize=font0)
blevels = [-4500, -4000, -3500, -3000, -2500, -2000, -1500, -1000,-500, 0] # define levels for plotting (transition of colorbar)
cbar.set_ticks(blevels)
cbar.set_ticklabels(blevels)
#grid
ax1.grid(linewidth=0.5)
ax2.grid(linewidth=0.5)
ax3.grid(linewidth=0.5)
# moving third figure slightly to the right
pos = ax3.get_position()
pos.x0 += 0.2
pos.x1 += 0.2
ax3.set_position(pos)
# moving second figure slightly to the left
pos = ax2.get_position()
pos.x0 -= 0.02
pos.x1 -= 0.02
ax2.set_position(pos)
# make plot boarders gray
for ax in [ax1, ax2, ax3]:
for spine in ax.spines.values():
spine.set_edgecolor('grey')
plt.tight_layout()
# save
fig.savefig('figures/MS1_scenarios.png', dpi=300, bbox_inches='tight')
fig.savefig('figures/MS1_scenarios.pdf', dpi=300, bbox_inches='tight')