-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMapsCoralCoverClean.m
339 lines (296 loc) · 12.7 KB
/
MapsCoralCoverClean.m
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
%% Make Maps
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evolutionary model for coral cover (from Baskett et al. 2009) %
% modified by Cheryl Logan ([email protected]) %
% last updated: 1-6-15 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = MapsCoralCoverClean(fullDir, Reefs_latlon, activeReefs, ...
lastYearAlive, events85_2010, eventsAllYears, frequentBleaching, ...
mortState, bleachState, fullYearRange, modelChoices, C_yearly, C_seed, con)
% Add paths and load mortality statistics
%load(strcat('~/Dropbox/Matlab/SymbiontGenetics/',filename,'/201616_testNF_1925reefs.mat'),'Mort_stats')
format shortg;
% filename = '201616_figs'; %filename = strcat(dateString,'_figs'); mkdir(filename); % location to save files
% map %% NOTE: worldmap doesn't seem to be working on work computer
%%
% Ten files are output, most of which are never used. Set this true to produce
% only the one map which is used in publication figures. This is the slowest
% part of postprocessing, so save time and disk space!
healthyV2Only = true;
%% yearRange provides the scale for all time related color bars. It spans
% all years during which a reef died, rounded out to the nearest ten. If
% no reefs die, a default of 1960 to 2100 is used.
if ~any(lastYearAlive)
yearRange = [1960 2100];
else
yearRange = [min(lastYearAlive) max(lastYearAlive)];
% Plotting chokes if the values are equal.
if yearRange(1) == yearRange(2)
yearRange(2) = yearRange(2) + 1;
end
% Also round out to a multiple of 10
if mod(yearRange(1), 10)
yearRange(1) = 10*floor(yearRange(1)/10);
end
if mod(yearRange(2), 10)
yearRange(2) = 10*ceil(yearRange(2)/10);
end
end
% We need to map a spot for all reefs, to show those that never bleached.
% Not every reef has a last mortality, but all have BLEACH8510 stats.
activeLatLon(1:length(activeReefs), 1) = Reefs_latlon(activeReefs, 1);
activeLatLon(1:length(activeReefs), 2) = Reefs_latlon(activeReefs, 2);
% A scale for "red=bad, blue=good" plots.
customColors = customScale();
if ~healthyV2Only
%% Make map of last mortality event recorded
tName = strcat(modelChoices,'. Year Corals No Longer Persist');
fileBase = strcat(fullDir, modelChoices,'_LastYrCoralMap');
% Green points everywhere
oneMap(12, activeLatLon(:, 1), activeLatLon(:, 2), [0 0.8 0], yearRange, parula, tName,'', false);
% Color-scaled points where there is a last year
outFile = strcat(fileBase, '.pdf');
if any(lastYearAlive)
ind = find(lastYearAlive);
oneMap(12, Reefs_latlon(ind, 1), Reefs_latlon(ind, 2), lastYearAlive(ind), yearRange, customColors, tName, outFile, true);
end
% This one may be post-processed, so save .fig
if verLessThan('matlab', '8.2')
saveas(gcf, fileBase, 'fig');
else
savefig(strcat(fileBase,'.fig'));
end
%% Make map showing # all bleaching events bn 1985-2010
tName = 'Bleaching Events Between 1985-2010';
fileBase = strcat(fullDir, modelChoices,'_MortEvents8510Map');
outFile = strcat(fileBase, '.pdf');
oneMap(13, activeLatLon(:, 1), activeLatLon(:, 2), events85_2010(activeReefs), [], jet, tName, outFile, false);
% Another one with postprocessing...
if verLessThan('matlab', '8.2')
saveas(gcf, fileBase, 'fig');
else
savefig(strcat(fileBase,'.fig'));
end
%% Figure 14 Make map showing # all bleaching events
rangeText = sprintf('%d-%d',fullYearRange);
tName = strcat('Bleaching Events Between ', rangeText);
outFile = strcat(fullDir, modelChoices,'_AllMortEventsMap','.pdf');
oneMap(14, activeLatLon(:, 1), activeLatLon(:, 2), eventsAllYears(activeReefs), [], jet, tName, outFile, false);
%% Figure 15 Same as 14 but with restricted color scale
cRange = [0, 20];
outFile = strcat(fullDir, modelChoices,'_AllMortEventsMap_Scale20','.pdf');
oneMap(15, activeLatLon(:, 1), activeLatLon(:, 2), eventsAllYears(activeReefs), cRange, jet, tName, outFile, false);
%% Figure 16 Same as 14 but with log2 of the number of events
% tName = 'Bleaching Events Between 1861-2100 (log base 2)';
% outFile = strcat(fullDir, modelChoices, '_AllMortEventsMap_log2', '.pdf');
% oneMap(16, activeLatLon(:, 1), activeLatLon(:, 2), log2(eventsAllYears(activeReefs)), [], jet, tName, outFile, false);
%% Figure 17. Maps first year of unhealthy coral, defined as:
% - no full-reef bleaching
% - no full-reef mortality
% - not currently bleached
% This can be expressed as the minimum the first year for each of those
% indicators.
% Store indexes, not years in lastHealthy, until just before plotting.
firstUnhealthy = NaN(length(Reefs_latlon), 1);
fbMass = frequentBleaching(:, :, 1);
msMass = mortState(:, :, 1);
bBoth = bleachState(:, :, end);
for k = activeReefs
% Frequent
ind = find(fbMass(k, :), 1, 'first');
if ~isempty(ind)
firstUnhealthy(k) = ind;
end
% Mortality (It may be that bleaching is always flagged when this is
% true, so it could be skipped - but for now be safe.)
ind = find(msMass(k, :, 1), 1, 'first');
if ~isempty(ind)
firstUnhealthy(k) = min(firstUnhealthy(k), ind);
end
% Current bleaching
ind = find(bBoth(k, :), 1, 'first');
if ~isempty(ind)
firstUnhealthy(k) = min(firstUnhealthy(k), ind);
end
end
% Convert from indices to year. NaN stays NaN.
firstUnhealthy = firstUnhealthy + fullYearRange(1) - 1;
tName = strcat(modelChoices,'. First Year of Unhealthy Reef');
fileBase = strcat(fullDir, modelChoices, '_FirstUnHealthyReef');
outFile = strcat(fileBase, '.pdf');
oneMap(17, activeLatLon(:, 1), activeLatLon(:, 2), firstUnhealthy(activeReefs), [], customColors, tName, outFile, false);
%% Figure 18. Maps last year of healthy coral, defined as:
% - no full-reef bleaching
% - no full-reef mortality
% - not currently bleached
% Do this by combining all the flags and then looking for the last year
% of health.
%
% Dimensions:
% frequentBleaching, mortState, and bleachState are all reefs x years x coral types.
% mortState and bleachState include an extra column for "all"
% Store indexes, not years in lastHealthy, until just before plotting.
lastHealthy = NaN(length(Reefs_latlon), 1);
combo = frequentBleaching(:, :, 1) | mortState(:, :, 1) | bleachState(:, :, 1);
% Now we need do find the last time the value is false (healthy)
for k = activeReefs
ind = find(~combo(k, :), 1, 'last');
if ~isempty(ind)
lastHealthy(k) = ind;
end
end
% Convert from indices to year. NaN stays NaN.
lastHealthy = lastHealthy + fullYearRange(1) - 1;
lastYearRange = [1950 2100];
tName = strcat(modelChoices,'. Last Year of Healthy Reef');
fileBase = strcat(fullDir, modelChoices, '_LastHealthyReef');
outFile = strcat(fileBase, '.pdf');
oneMap(18, activeLatLon(:, 1), activeLatLon(:, 2), lastHealthy(activeReefs), lastYearRange, customColors, tName, outFile, false);
% This one may be post-processed, so save .fig
if verLessThan('matlab', '8.2')
saveas(gcf, fileBase, 'fig');
else
savefig(strcat(fileBase,'.fig'));
end
end
%% Figure 20. Maps last year of healthy coral, defined as:
%
% Fig 19 showed early problems when massive corals were still healthy. Try
% basing more on just the massives.
%
% In terms of the events flagged:
% - bleaching - no longer included?
% - no mortality of either coral type
% - no high-frequency bleaching of either coral type
% Do this by combining all the flags and then looking for the last year
% of health.
%
% Dimensions:
% frequentBleaching, mortState, and bleachState are all reefs x years x coral types.
% mortState and bleachState include an extra column for "all"
% Store indexes, not years in lastHealthy, until just before plotting.
lastHealthy = NaN(length(Reefs_latlon), 1);
% & operator has precedence before |
combo = frequentBleaching(:, :, 1) | (mortState(:, :, 1) & mortState(:, :, 2));
% Now we need to find the last time the value is false (healthy)
for k = activeReefs
ind = find(~combo(k, :), 1, 'last');
if ~isempty(ind)
lastHealthy(k) = ind;
end
end
% Convert from indices to year. NaN stays NaN.
lastHealthy = lastHealthy + fullYearRange(1) - 1;
lastYearRange = [1950 2100];
tName = strcat(modelChoices,'. Last Year of Healthy Reef');
fileBase = strcat(fullDir, modelChoices, '_LastHealthyBothTypesV2');
outFile = strcat(fileBase, '.pdf');
oneMap(20, activeLatLon(:, 1), activeLatLon(:, 2), lastHealthy(activeReefs), lastYearRange, customColors, tName, outFile, false);
% This one may be post-processed, so save .fig
if verLessThan('matlab', '8.2')
saveas(gcf, fileBase, 'fig');
else
savefig(strcat(fileBase,'.fig'));
end
if healthyV2Only
return;
end
%% Figure 21. Maps last year of healthy coral, defined as:
%
% This is similar to Figure 20, but requiring corals to be at 4 * Seed for at
% least one coral type to be considered healthy.
%
% Dimensions:
% frequentBleaching, mortState, and bleachState are all reefs x years x coral types.
% mortState and bleachState include an extra column for "all"
% Store indexes, not years in lastHealthy, until just before plotting.
lastHealthy = NaN(length(Reefs_latlon), 1);
% There's a catch when not all reefs are computed. Most arrays here have the
% full 1925 reef rows or columns, but C_yearly is only for the computed reefs.
% We'll need an expanded version in that case.
if size(C_yearly, 2) < size(Reefs_latlon, 1)
% C_yearly is years/reefs/corals, but the other arrays are
% reefs/years/corals!
C_full = zeros(size(frequentBleaching));
num = 0;
for k = activeReefs
num = num + 1;
C_full(k, :, :) = C_yearly(:, num, :);
end
else
C_full = permute(C_yearly, [2 1 3]);
end
adequateCover = false(size(C_full));
%adequateCover(:, :, 1) = C_full(:, :, 1) > C_seed(1) * 4;
%adequateCover(:, :, 2) = C_full(:, :, 2) > C_seed(2) * 4;
adequateCover(:, :, 1) = C_full(:, :, 1) > 0.2 * con.KCm;
adequateCover(:, :, 2) = C_full(:, :, 2) > 0.2 * con.KCm;
% Compute combo. 1 is bad.
% & operator has precedence before |
combo = (~adequateCover(:, :, 1) & ~adequateCover(:, :, 2)) | (frequentBleaching(:, :, 1) | (mortState(:, :, 1) & mortState(:, :, 2)));
% Now we need to find the last time the value is false (healthy)
for k = activeReefs
ind = find(~combo(k, :), 1, 'last');
if ~isempty(ind)
lastHealthy(k) = ind;
end
end
% Convert from indices to year. NaN stays NaN.
lastHealthy = lastHealthy + fullYearRange(1) - 1;
lastYearRange = [1950 2100];
tName = strcat(modelChoices,'. Last Year of Healthy Reef');
fileBase = strcat(fullDir, modelChoices, '_LastHealthyBothTypesV2');
outFile = strcat(fileBase, '.pdf');
oneMap(21, activeLatLon(:, 1), activeLatLon(:, 2), lastHealthy(activeReefs), lastYearRange, customColors, tName, outFile, false);
% This one may be post-processed, so save .fig
if verLessThan('matlab', '8.2')
saveas(gcf, fileBase, 'fig');
else
savefig(strcat(fileBase,'.fig'));
end
end % End the main MapsCoralCover function.
% Arguments:
% n figure number
% lons longitudes (was [events.lon])
% lats latitudes
% values what to plot at each position
% cRange man and max data values for color scale
% t title
% outFile pdf output file
function [] = oneMap(n, lons, lats, values, cRange, cMap, t, outFile, add)
f = figure(n);
if add
hold on;
else
clf;
% first pass only:
%m_proj('miller', 'lat', [-40 40]); % , 'lon', 155.0); - offsets map, but drops some data!
% Center on the coral triangle, skipping coral-free longitudes around
% zero.
m_proj('miller', 'lat',[-40 40],'long',[20 340]); % [0 360] for world, but no reefs from -28.5 to +32 longitude
m_coast('patch',[0.7 0.7 0.7],'edgecolor','none');
m_grid('box','off','linestyle','none','backcolor',[.9 .99 1], ...
'xticklabels', [], 'yticklabels', [], 'ytick', 0, 'xtick', 0);
end
% Shift since maps is now in positive longitude terms.
idx = find(lons < 0);
lons(idx) = lons(idx) + 360; % for shifted map (0 to 360 rather than -180 to 180)
% Points with last-year mortality values:
[LONG,LAT] = m_ll2xy(lons,lats); hold on % convert reef points to M-Map lat long
scatter(LONG,LAT,5, values) ; % plot bleaching events onto map
if isempty(cMap)
colormap default;
else
colormap(cMap)
end
if ~isempty(cRange)
caxis(cRange);
end
colorbar
title(t)
if ~isempty(outFile)
print(f, '-dpdf', '-r200', outFile);
end
hold off;
end