close all;
clear all;
disp(' ')
disp(' ')
disp(' John Horel')
disp(' ATMOS 5040/ Spring 2017')
disp(' Chapter 2f')
disp('requires loadjson.m')
disp('requires varargin2struct.m')
disp('jsonopt.m')

%get time series of wind speed and direction from WBB for a month

%documentation on MesoWest api
%http://synopticlabs.org
%June 2016
api_string = strcat('http://api.synopticlabs.org/v2/stations/timeseries?stid=wbb&start=201606010000&end=201606302355&vars=wind_speed,wind_direction,wind_gust&obtimezone=local&token=tk2017mar0413y8nZBC35uyjSYMJqVyYfY')
% January 2017
%api_string = strcat('http://api.synopticlabs.org/v2/stations/timeseries?stid=wbb&start=201701010000&end=201701312355&vars=wind_speed,wind_direction,wind_gust&obtimezone=local&token=tk2017mar0413y8nZBC35uyjSYMJqVyYfY')
%download json result as a string
api_return = urlread(api_string);
% interpret the json string using the loadjson function
json = loadjson(api_return);
% json is a structure that contains UNITS, STATION, and SUMMARY
% assign units of the variables. 
units = json.UNITS;
% obtain summary associated with the api call
summary = json.SUMMARY;
% parse the STATION info. STATION is a cell vector inside of the json structure
% push cell vector into a structure: station
station = cell2mat(json.STATION);
%get all the times, strip out the T and Z use string replace, save as vector and matlab time
times = station.OBSERVATIONS.date_time;
time = strrep(times,'T',' ');
tim = strrep(time,'Z','');
tv=datevec(tim,'yyyy-mm-dd HH:MM:SS');
td=datenum(tim,'yyyy-mm-dd HH:MM:SS');

% get the speed and direction data
%set_1 refers to that a site could have more 
%than one sensor reporting the same variable
spd = station.OBSERVATIONS.wind_speed_set_1;
gst = station.OBSERVATIONS.wind_gust_set_1;
dirc = station.OBSERVATIONS.wind_direction_set_1;
% there are some missing directions, which is awkward
% so have to search and in this case replace with 0
dire=find(cellfun(@isempty,dirc));
%brute force replace empty with 0. Not a recommended solution
for i=1:length(dire)
    dirc{1,dire(i)}=0;
end
dir= cell2mat(dirc);

%convert spd,dir to u,v. matlab- 0 angle is relative to x axis
% and counterclockwise
u = spd .* sind(dir - 180);
v = spd .* cosd(dir - 180);

figure(1)
subplot(2,2,1)
sp1= 0:0.5:15;
hist(spd,sp1);
axis([0,15,0,1000]);
hold on;
xlabel('Wind Speed (m/s)');
ylabel('Count');
title('Histogram: John Horel 02/02/2017');

subplot(2,2,2)
cdfplot(spd)
xlabel('Wind Speed (m/s)');
ylabel('Cumulative Distributions');
title('Cumulative Distributions: John Horel 02/02/17');

subplot(2,2,3)
dir1= 0:5:360;
hist(dir,dir1);
axis([0,360,0,400]);
xlabel('Wind Direction');
ylabel('Count');


subplot(2,2,4)
cdfplot(dir)
xlabel('Horizontal Wind Direction');
ylabel('Cumulative Distribution');

%plot a 2 panel plot of
%top- histogram of wind gust from 0-15 at intervals of 0.5 m/s
%bottom- cdf of wind gust


%compute mean, median, mode, sample std, unbiased estimate of the population standard deviation, and skewness of gust


figure(3)
u1= -5:0.25:5;
subplot(2,2,1)
hist(u);
hold on;
xlabel('Zonal Wind Speed (m/s)');
ylabel('Count');
title('Histograms: John Horel 02/02/17');

subplot(2,2,2)
cdfplot(u);
xlabel('Zonal Wind (m/s)');
ylabel('Cumulative Distribution');
title('Cumulative Distribution: John Horel 02/02/2017');

subplot(2,2,3)
v1= -10:0.5:10;
hist(v,v1);
hold on;
xlabel('Meridional Wind (m/s)');
ylabel('Count');

subplot(2,2,4)
cdfplot(v);
xlabel('Meridional Wind(m/s)');
ylabel('Cumulative Distribution');

figure(4)
% wind rose
n=length(dir);
ctrl = 0;
ctrm = 0;
%loop over all elements
for i=1:n
 %convert direction to radians  
 dirr(i,1) = dir(i) * pi/180;
 %count light winds first
    if (spd(i) < 2)
        ctrl=ctrl+1;
    dirl(ctrl,1) = dirr(i,1);
    end
    %count medium winds next
if (spd(i)<4)
    ctrm=ctrm+1;
    dirm(ctrm,1) = dirr(i,1);
end
end
%30 degree increments
x = [pi/12:pi/6:2*pi-pi/12];
%plot all winds first to set axes
h=rose(dirr,x);
set(h,'LineWidth',4,'Color','k');
hold on;
%plot medium winds
hm=rose(dirm,x);
set(hm,'LineWidth',2,'Color','r');
%plot light winds
hl=rose(dirl,x);
set(hl,'LineWidth',3,'Color','b');
title('WBB Hodograph: John Horel 02/02/17');
%by default the rose is reversed relative to 
%convention used in the atmospheric sciences
set(gca,'View',[-90 90],'YDir','reverse');

ctr = zeros(1,24);
%loop over all elements of time series
for i=1:n
 %load up arrays for each hour of the day
 h = tv(i,4)+1;
 ctr(1,h)=ctr(1,h)+1;
 sph(ctr(1,h),h)=spd(i);
 dh(ctr(1,h),h)=dir(i);
 uh(ctr(1,h),h)=u(i);
 vh(ctr(1,h),h)=v(i);
end
%compute mean wind speed for each hour
spm = sum(sph)./ctr;
% compute mean components each hour
um = sum(uh)./ctr;
vm = sum(vh)./ctr;
% compute resultant wind speed for each hour
ressp = sqrt(um.*um+vm.*vm);
hrlbl = 0:1:23;
figure(5)
hold on
plot(hrlbl,spm,'r')
plot(hrlbl,ressp,'g')
axis([0 23 0 5])
xlabel('Hour of day');
ylabel('Wind Speed (m/s)');
title('WBB Mean and Resultant Wind Speed: John Horel 02/02/2017');
