-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSolana_RadiationStress_MainScript.m
118 lines (76 loc) · 2.43 KB
/
Solana_RadiationStress_MainScript.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
%% Solana Radiaiton Stress Investigation
% Author: Trenton Saunders
% Date: 06-28-2024
% Calculate alongshore radiation stress gradient (Ludka et al. 2018)
close all
clear all
clc
dir_string = cd;
%% Set Up NCTOOLBOX
cd ..\
addpath(fullfile(cd,'nctoolbox'))
setup_nctoolbox
cd(dir_string)
%% USER ENTERS STATION NUMBER AND START/END DATES FOR PLOT
MOPstn = 675:-1:630;
% startdate = '01/01/2021 00:00';
% enddate = '01/31/2021 23:59';
%% Pull MOP data from THREDDS
for i = 1:length(MOPstn)
test_url1 = 'https://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/model/MOP_alongshore/D0';
test_url2 = '_hindcast.nc';
test_url = [test_url1,num2str(MOPstn(i)),test_url2];
MOP_info = ncinfo(test_url)
Data.Sxy = ncread(test_url,'waveSxy');
Data.Hs = ncread(test_url,'waveHs');
Data.Time = ncread(test_url,'waveTime');
Data.MOP(i) = MOPstn(i);
Data.Lat(i) = ncread(test_url,'metaLatitude');
Data.Lon(i) = ncread(test_url,'metaLongitude');
Data.Average_Sxy(i) = mean(Data.Sxy);
Data.Average_Hs(i) = mean(Data.Hs);
[Data.X(i),Data.Y(i),Data.utmzone{i}] = deg2utm(Data.Lat(i),Data.Lon(i));
end
%%
% figure
% plot(Data.MOP(:),Data.Average_Sxy(:))
%% Figure
figure
subplot(1,2,1)
geobubble(Data.Lat(:),Data.Lon(:),Data.Average_Hs(:))
geobasemap satellite
sgtitle('$\overline {H_{s}}$ ', 'Interpreter','Latex')
subplot(1,2,2)
geobubble(Data.Lat(:),Data.Lon(:),Data.Average_Sxy(:))
geobasemap satellite
sgtitle('$\overline {S_{xy}}$ ', 'Interpreter','Latex')
%% Calculate mean Sxy alongshore gradient (nned to check this calculation)
figure
plot(Data.Y,Data.Average_Sxy,'.')
diff_Y = diff(Data.Y)
diff_Sxy = diff(Data.Average_Sxy)
grad_Sxy = diff(Data.Average_Sxy)./diff(Data.Y);
%% Calculate corresponding location of gradinet calcultion
meanX = movmean(Data.X,2,"Endpoints","discard");
meanY = movmean(Data.Y,2,"Endpoints","discard");
Sxy_Range = linspace(min(grad_Sxy)*10000,max(grad_Sxy)*10000,100);
Sxy_Range = linspace(-.4,.4,100);
colors = parula(100)
figure
for j = 1:length(meanX)
[MeanLat(j),MeanLon(j)] = utm2deg( meanX(j),meanY(j),'11 S' );
[~, indx] = min(abs(Sxy_Range - 10000*grad_Sxy(j)))
geoscatter(MeanLat(j),MeanLon(j),200,[colors(indx,:)],'filled')
hold on
end
geobasemap satellite
colorbar
clim([Sxy_Range(1),Sxy_Range(end)])
%%
Sxy_Range = linspace(min(grad_Sxy)*10000,max(grad_Sxy)*10000,100);
cool(100)
% Data.DateTime = datetime(Data.Time,'ConvertFrom','posixtime','TimeZone','UTC');
%
% figure
% plot(Data.DateTime,Data.Sxy)
%