function [X,Y,AVGSEC,BASELINE,NUMPOINTS,MEANTEMP,STDTEMP,MEANSAL,STDSAL,MEANST,STDST,MEANC,STDC]=cl_climo_xshelf(filename,season,bathyline,bathychoice,maxdepth,onshore,offshore,hbinsize,vbinsize,westbound,eastbound,southbound,northbound,minnumpts,qc,decim,smoothing,plotyes); % % CL_CLIMO_XSHELF This program takes a file of unevenly spaced hydrographic % data and creates cross-shelf climatological mean fields based on user-specified % inputs. Users can also elect to quality-control and decimate the data. % % function [X,Y,AVGSEC,BASELINE,NUMPOINTS,MEANTEMP,STDTEMP,MEANSAL,STDSAL,MEANST,STDST,MEANC,STDC]=cl_climo_xshelf(filename,season,bathyline,bathychoice,maxdepth,onshore,offshore,hbinsize,vbinsize,westbound,eastbound,southbound,northbound,minnumpts,qc,decim,smoothing,plotyes); % % Supporting m-files: % - The m_map toolbox, free at this website: http://www2.ocgy.ubc.ca/~rich/map.html % - If you would like smoothed and extrapolated output, you will need to % download Roger Goldsmith's adaptation of the PlotPlus ppzgrid software % (all Matlab m-files). Smoothing is optional--if you don't want smoothed % output you don't need to download this toolbox. % % Inputs: % filename - specifies name of data file % Data file, in ASCII text format, must have the following columns: % 1. Decimal longitude % 2. Decimal latitude % 3. Year % 4. Month % 5. Depth (meters) % 6. Temperature (C) % 7. Salinity (PSU) % season string - is a listing of all seasons (by number) to include in average. For example, Jan & Feb would be '1 2'. Enter 'All' for all. % bathyline BASELINE isobath for the climatology in meters. % bathychoice either 0, which forces the user to hand-pick the BASELINE, or a 2 column array of lon/lat pairs defining the BASELINE, which bypasses the BASELINE chooser routine. % maxdepth maximum depth of the climatology in meters. % onshore extent of onshore domain in kilometers. % offshore extent of offshore domain in kilometers. % hbinsize horizontal bin size in kilometers. % vbinsize vertical bin size in meters. % westbound western domain bound in decimal degrees. % eastbound eastern domain bound in decimal degrees. % southbound southern domain bound in decimal degrees. % northbound northern domain bound in decimal degrees. % minnumpts minimum number of data points to use for an average. If this number is not met, a NaN will result for the mean for that grid cell. % qc 0 or 1; 0 will skip the quality control routine and 1 will run it. % decim 0 or 1; 0 will bypass decimation routine and 1 will run it. % smoothing smoothing applied by ppzgrid routine. A value if 2 is typical for hydrographic section data. Enter 0 to bypass smoothing. % plotyes 0, 1, or 2; 0 will suppress plots, 1 will show plots, and 2 will show plots and save them out as PNG files to the current directory. % % Outputs: % X Cross-shelf distance of output grid points in KM (meshgrid format). % Y Depth of output grid points in M (meshgrid format). % AVGSEC 2-column file defining the average cross-shelf bathymetry profile. The first column is distance in KM and the second is depth in M. % BASELINE 2-column file of the user-selected BASELINE. This is preserved as an output variable so that for future runs, you can use this variable as an INPUT?see bathychoice, above. % NUMPOINTS Number of points per output grid node. This, like all of the other output fields, are in the meshgrid format for easy plotting. For example, pcolor(X,Y,NUMPOINTS). % MEANTEMP Mean temperature (C). % STDTEMP Standard deviation of temperature. % MEANSAL Mean salinity (PSU). % STDSAL Standard deviation of salinity. % MEANST Mean sigma-theta (density) (kg/m3) (sigma-theta and sound speed outputs computed using a 1024-realization monte-carlo routine with mean & std temp/salinity as inputs). % STDST Standard deviation of sigma-theta. % MEANC Mean sound speed (m/s). % STDC Standard deviation of sound speed. % % Example: Nantucket Shoals test case % [X,Y,AVGSEC,BASELINE,NUMPOINTS,MEANTEMP,STDTEMP,MEANSAL,STDSAL,MEANST,STDST,MEANC,STDC]=cl_climo_xshelf('ns-all-data.dat','8 9',100,0,200,60,40,10,10,-72,-69,39,41,10,1,1,2,1); % % Given a data file 'ns-all-data.dat', a climatology will be created for August % and September data, along a BASELINE isobath of 100m (user must choose basline % using mouse), extending from 60km onshore to 40km offshore and 200m deep. % The vertical bin size is 10 meters and the horizontal bin size is 10km. The % study area bounds are 72-69 West, 39-41 North. A minimum of 10 realizations % must be present to make a mean, the QC and decimation routines will be called, % a factor of 2 smoothing will be applied to the results, and the data will % be plotted but the plots will not be saved out. % % Version 4.0, 11Oct2006, Chris Linder clinder@whoi.edu % http://www.whoi.edu/science/PO/people/clinder/software_climo_xshelf.html % % This program is free software; you can redistribute it and/or % modify it under the terms of the GNU General Public License % as published by the Free Software Foundation; either version 2 % of the License, or (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. %%%%%%%%%%%%%%%%%%%%%%%%%% % Create study area plot % %%%%%%%%%%%%%%%%%%%%%%%%%% warning off MATLAB:divideByZero more off; if plotyes==2 saveprefix=[filename(1:end-4),'_']; end; disp(' '); eval(sprintf('disp(''Loading %s...'');',filename)); pause(0.1); eval(sprintf('rawdata=load(''%s'');',filename)); % Extract 5 min Terrainbase bathymetry file for this area (m_map toolbox) m_proj('mercator','lon',[westbound eastbound],'lat',[southbound northbound]); [data,LON,LAT]=m_tbase([westbound-1 eastbound+1 southbound-1 northbound+1]); % Extract data for chosen season if season(1)=='A' | season(1)=='a' numseasons=[1:1:12]; else numseasons=str2num(season); end; i1=find(rawdata(:,4)==1); i2=find(rawdata(:,4)==2); i3=find(rawdata(:,4)==3); i4=find(rawdata(:,4)==4); i5=find(rawdata(:,4)==5); i6=find(rawdata(:,4)==6); i7=find(rawdata(:,4)==7); i8=find(rawdata(:,4)==8); i9=find(rawdata(:,4)==9); i10=find(rawdata(:,4)==10); i11=find(rawdata(:,4)==11); i12=find(rawdata(:,4)==12); ichoice=[]; for sealoop=1:size(numseasons,2) eval(sprintf('ichoice=[ichoice; i%d];',numseasons(sealoop))); end; seadist=[]; ALLPOINTS=rawdata(ichoice,:); % Convert depth column to positive numbers ALLPOINTS=[ALLPOINTS(:,1:4) abs(ALLPOINTS(:,5)) ALLPOINTS(:,6:7)]; BASELINE=bathychoice; if plotyes | length(bathychoice)==1 f_map=figure; % Grid setup m_proj('mercator','lon',[westbound eastbound],'lat',[southbound northbound]); m_gshhs_i('patch',[.7 .7 .7]); m_grid('box','fancy'); hold on; m_plot(ALLPOINTS(:,1),ALLPOINTS(:,2),'g.'); [cb,hb]=m_tbase('contour',[-bathyline -bathyline],'b'); hold on; set(hb,'linewidth',2,'linecolor','b'); [cball,hball]=m_tbase('contour',[-4000 -2000 -1000 -500 -200 -100 -50],'k'); hold on; for fixloop=1:length(hball) set(hball(fixloop),'linecolor','k'); end; clabel(cb,hb,'color','black','fontsize',8); clabel(cball,hball,'color','black','fontsize',8); title('Choose points with LMB. Click RMB when done.'); if(length(bathychoice)==1) disp(' '); disp('This plot shows your study area with data points in green.'); disp('The bathymetry is contoured in black. Your initially chosen bathymetry'); disp('BASELINE is shown as a thicker blue line. You must now select the BASELINE'); disp('for your climatology using the mouse. Select points along your chosen '); disp('bathymetry line (shown in bold, blue). As you click, the points will be'); disp('plotted on the map in red. Click the right mouse button when done.'); BASELINE=[]; [x,y,button]=ginput(1); while button ~= 3 [LONGITUDE,LATITUDE]=m_xy2ll(x,y); BASELINE=[BASELINE; LONGITUDE, LATITUDE]; m_plot(BASELINE(:,1),BASELINE(:,2),'r.'); m_plot(BASELINE(:,1),BASELINE(:,2),'r-','linewidth',2); [x,y,button]=ginput(1); end; end; end; % Now we have a nice user-defined BASELINE. Need to figure out which is "onshore" and which is "offshore" % First, check the depth of a point 30 degrees to the right and half of the total BASELINE distance % from the first point chosen... BASELINE_azimuth=azimuth(BASELINE(1,2), BASELINE(1,1), BASELINE(end,2), BASELINE(end,1)); BASELINE_distance=deg2km(distance(BASELINE(1,2), BASELINE(1,1), BASELINE(end,2), BASELINE(end,1))); [checklat,checklon]=reckon(BASELINE(1,2),BASELINE(1,1),km2deg(BASELINE_distance/2),BASELINE_azimuth+30); checkelev=interp2(LON,LAT,data,checklon,checklat); % if checkelev is deeper than bathyline, then the offshore direction is to % the right of the first BASELINE point looking down the line... % Build study box... onshore=abs(onshore); offshore=abs(offshore); if checkelev < -bathyline d1=offshore; d2=onshore; else d1=onshore; d2=offshore; end; studybox=[]; profiles=[]; % Offshore polygon bound first for boxloop=1:size(BASELINE,1) % start in a direction ANTICLOCKWISE from BASELINE start [pt1lat,pt1lon]=reckon(BASELINE(boxloop,2),BASELINE(boxloop,1),km2deg(d1),BASELINE_azimuth+90); studybox=[studybox; pt1lon pt1lat]; end; studybox_offshore=[studybox; flipud(BASELINE)]; % Onshore polygon next studybox_onshore=BASELINE; for boxloop=1:size(BASELINE,1) % add points continuing AROUND CCW (for m_patch command) index=size(BASELINE,1)-boxloop+1; [pt2lat,pt2lon]=reckon(BASELINE(index,2),BASELINE(index,1),km2deg(d2),BASELINE_azimuth-90); profile=[]; for profileloop=1:21 dstep=km2deg((d1+d2)/20 * (profileloop-1)); [p1,p2]=reckon(pt2lat,pt2lon,dstep,BASELINE_azimuth+90); profile=[profile; interp2(LON,LAT,data,p2,p1)]; end; studybox=[studybox; pt2lon pt2lat]; studybox_onshore=[studybox_onshore; pt2lon pt2lat]; profiles=[profiles profile]; end; studybox=[studybox; studybox(1,:)]; profile_x=-onshore:(d1+d2)/20:offshore; if checkelev > -bathyline profiles=flipud(profiles); end; if plotyes > 0 | length(bathychoice)==1 m_patch(studybox(:,1),studybox(:,2),[1 .7 .7]); m_grid('box','fancy'); set(findobj('tag','m_grid_color'),'facecolor','none'); title(['Months: ',season,'; RED=study area, BLK=bathymetry, BLUE=BASELINE isobath, GRN=data']); % Replot part of the map so that the study area appears underneath m_plot(ALLPOINTS(:,1),ALLPOINTS(:,2),'g.'); [cb,hb]=m_tbase('contour',[-bathyline -bathyline],'b'); hold on; set(hb,'linewidth',2,'linecolor','b'); [cball,hball]=m_tbase('contour',[-4000 -2000 -1000 -500 -200 -100 -50],'k'); hold on; for fixloop=1:length(hball) set(hball(fixloop),'linecolor','k'); end; clabel(cb,hb,'color','black','fontsize',8); clabel(cball,hball,'color','black','fontsize',8); m_plot(BASELINE(:,1),BASELINE(:,2),'r.'); m_plot(BASELINE(:,1),BASELINE(:,2),'r-','linewidth',2); end; % Create an average bathymetry section along BASELINE. % First, use reckon to get 20 points in a perpendicular line across each % user-selected point. Use interp2 to get the depth from Terrainbase at % each point. Then, add each profile into an array. Finally, % average the entire array into 1 profile. Rename avgsec. disp(' '); disp('Creating average bathymetry section...'); AVGSEC=[profile_x' nanmean(profiles,2)]; if plotyes > 0 f_bathy=figure; plot(profile_x,profiles,'k'); xlim([-onshore offshore]); ylim([-maxdepth 0]); hold on; plot(AVGSEC(:,1),AVGSEC(:,2),'r','linewidth',2); xlabel('Distance (km)'); ylabel('Depth (m)'); yl=get(gca,'ytick'); set(gca,'ytick',yl); set(gca,'yticklabel',abs(yl)); title('Cross-shelf bathymetry realizations (black) and mean (red)'); if plotyes==2 print('-dpng',[saveprefix 'bathy.png']); end; end; % Go through data file, extracting casts inside the study area and using % interp1 to put each interpolated cast onto the regular (vbin) depth scale. disp(' '); disp('Extracting casts, computing derivative fields, interpolating onto depth bins, and sorting into cross-shelf bins... '); bins_dist=[-onshore + hbinsize/2 : hbinsize : offshore - hbinsize/2]; bins_depth=[-abs(vbinsize/2):-abs(vbinsize):-(abs(maxdepth)-abs(vbinsize/2))]; placeholder=1; allcasts=[]; for rownum=2:size(ALLPOINTS,1) % 1 possibility - lat OR lon OR month OR year changes OR depth decreases: % then it's a new cast... if ALLPOINTS(rownum,1) ~= ALLPOINTS(rownum-1,1) | ALLPOINTS(rownum,2) ~= ALLPOINTS(rownum-1,2) | ALLPOINTS(rownum,3) ~= ALLPOINTS(rownum-1,3) | ALLPOINTS(rownum,4) ~= ALLPOINTS(rownum-1,4) | ALLPOINTS(rownum,5) < ALLPOINTS(rownum-1,5) castprops=ALLPOINTS(placeholder:rownum-1,:); % Need to check for duplicates in the depth record (rare) and % eliminate them [b,i,j]=unique(castprops(:,5)); castprops=castprops(i,:); % Calculate distance from BASELINE dist_from_BASELINE=(deg2km(distance(BASELINE(:,2),BASELINE(:,1),ones(size(BASELINE(:,2)))*castprops(1,2),ones(size(BASELINE(:,2)))*castprops(1,1)))); [mindist,i_mindist]=min(dist_from_BASELINE); % Now we have the minimum (perpendicular) distance to the BASELINE % Which side (onshore/offshore) is it on? Use inside function. if inside(castprops(1,1),castprops(1,2),studybox_onshore(:,1),studybox_onshore(:,2)) mindist=-1*mindist; elseif inside(castprops(1,1),castprops(1,2),studybox_offshore(:,1),studybox_offshore(:,2)) mindist=mindist; else mindist=Inf; end; if mindist >= -onshore & mindist <= offshore & size(castprops,1) > 1 % Cast must be within bounds and have more than one point! profile_temp=interp1(-1*abs(castprops(:,5)),castprops(:,6),bins_depth)'; profile_sal=interp1(-1*abs(castprops(:,5)),castprops(:,7),bins_depth)'; profile_pressure=sw_pres(abs(castprops(:,5)),ones(size(castprops(:,5)))*castprops(1,2)); profile_density=interp1(-1*abs(castprops(:,5)),sw_dens(castprops(:,7),castprops(:,6),profile_pressure),bins_depth)'; profile_soundspeed=interp1(-1*abs(castprops(:,5)),sw_svel(castprops(:,7),castprops(:,6),profile_pressure),bins_depth)'; HBin=floor(mindist/hbinsize) + 7; numdate=datenum(castprops(1,3),castprops(1,4),15); regcast=[HBin; castprops(1,1); castprops(1,2); castprops(1,3); castprops(1,4); numdate; profile_temp; profile_sal; profile_density; profile_soundspeed]; allcasts=[allcasts regcast]; end; placeholder=rownum; end; end; % Remove exact duplicate casts... totalextracted=size(allcasts,2); allcasts=unique(allcasts','rows')'; totalremaining=size(allcasts,2); num_elim_dupes=totalextracted-totalremaining; disp(' '); disp([num2str(totalextracted),' casts extracted. ',num2str(num_elim_dupes),' duplicate casts removed.']); disp(['Total number of casts remaining: ',num2str(totalremaining),'.']); if(qc) % Quality control f_qc=figure; % Show the casts graphically disp(' '); disp('You have chosen to quality control the data. This will be accomplished by an iterative process. You will be shown all of cast T/S profiles, plus the mean and standard deviation envelope.'); disp(' '); disp('Next, you will choose a multiple of the standard deviation envelope -- casts within this envelope will be kept, and those outside will be flagged. If the T OR S is outside at any depth point in the profile, the entire cast will be flagged.'); disp(' '); disp('The flagged casts will be plotted in red. You will now be given the opportunity to accept this decision. If you are not satisfied, you will start the process over. After this initial flag decision, you will be asked to review each cast for inclusion/exclusion.'); originalcasts=allcasts; notsatisfied=1; while(notsatisfied) figure(f_qc); clf; subplot(121); hold on; plot(allcasts(7:6+size(bins_depth,2),:),bins_depth,'b.-'); title('Temperature casts'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); ylabel('Depth (m)'); xlabel('Temperature (C)'); subplot(122); hold on; plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),:),bins_depth,'b.-'); title('Salinity casts'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); ylabel('Depth (m)'); xlabel('Salinity'); % Compute mean profile and standard deviation and overlay on profile plots temp_mean=nanmean(allcasts(7:6+size(bins_depth,2),:),2); temp_std=nanstd(allcasts(7:6+size(bins_depth,2),:),0,2); figure(f_qc); subplot(121); plot(temp_mean,bins_depth,'k','linewidth',3); plot(temp_mean-temp_std,bins_depth,'k','linewidth',2); plot(temp_mean+temp_std,bins_depth,'k','linewidth',2); sal_mean=nanmean(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),:),2); sal_std=nanstd(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),:),0,2); figure(f_qc); subplot(122); plot(sal_mean,bins_depth,'k','linewidth',2); plot(sal_mean-sal_std,bins_depth,'k','linewidth',2); plot(sal_mean+sal_std,bins_depth,'k','linewidth',2); disp(' '); stddevnum=input('Choose multiple of standard deviation (does not have to be integer): '); if(isempty(stddevnum)) disp('Default number of "1" entered...'); stddevnum=1; end; t_envelope_minus=temp_mean-temp_std*stddevnum; t_envelope_plus=temp_mean+temp_std*stddevnum; s_envelope_minus=sal_mean-sal_std*stddevnum; s_envelope_plus=sal_mean+sal_std*stddevnum; disp(' '); disp('Plotting your desired envelope in magenta...'); subplot(121); plot(t_envelope_minus,bins_depth,'m','linewidth',2); plot(t_envelope_plus,bins_depth,'m','linewidth',2); subplot(122); plot(s_envelope_minus,bins_depth,'m','linewidth',2); plot(s_envelope_plus,bins_depth,'m','linewidth',2); % Now, flag potentially bad casts - ones in which any point is out % of line. To do this we will need to loop through all of the depths and flag % bad indices alltempdata=allcasts(7:6+size(bins_depth,2),:); allsaldata=allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),:); i_out_t=[]; i_out_s=[]; i_good=[]; for subloop=1:size(alltempdata,1) i_out_t=[i_out_t find(alltempdata(subloop,:) <= t_envelope_minus(subloop) | alltempdata(subloop,:) >= t_envelope_plus(subloop) )]; i_out_s=[i_out_s find(allsaldata(subloop,:) <= s_envelope_minus(subloop) | allsaldata(subloop,:) >= s_envelope_plus(subloop) )]; end; % Need to do a unique to get rid of multiples i_out=unique([i_out_t i_out_s]); if isempty(i_out)==0 % Some flagged subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out),bins_depth,'r.-'); title('Temperature (red=mean, green=std)'); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out),bins_depth,'r.-'); title('Salinity (red=mean, green=std)'); disp(' '); disp(['Your chosen envelope has found ',num2str(length(i_out)),' out of range casts out of ',num2str(size(allcasts,2)),' total casts.']); else % None flagged disp(' '); disp(['Your chosen envelope has found ',num2str(length(i_out)),' out of range casts out of ',num2str(size(allcasts,2)),' total casts.']); end; disp(' '); notsatisfied=input('Red (bad) casts are now plotted. Would you like to rerun (1), or continue (0, or just press enter): '); end; % while notsatisfied % Next, display casts one by one and ask if they should be included or % excluded if isempty(i_out)==0 % Some flagged disp(' '); i_out_final=[]; figure(f_qc); clf; subplot(121); hold on; plot(allcasts(7:6+size(bins_depth,2),:),bins_depth,'b.-'); plot(t_envelope_minus,bins_depth,'k','linewidth',2); plot(t_envelope_plus,bins_depth,'k','linewidth',2); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); ylabel('Depth (m)'); xlabel('Temperature (C)'); title('Temperature (blue=keep, red=reject)'); subplot(122); hold on; plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),:),bins_depth,'b.-'); title('Salinity (blue=keep, red=reject)'); plot(s_envelope_minus,bins_depth,'k','linewidth',2); plot(s_envelope_plus,bins_depth,'k','linewidth',2); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); ylabel('Depth (m)'); xlabel('Salinity'); % Plot the first questionable cast outloop=1; subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out(outloop)),bins_depth,'m-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out(outloop)),bins_depth,'m-','linewidth',2); if size(i_out,2) > 1 % More than 1 questionable cast keepit=input('Kill all questionable casts (2), keep the magenta cast (1), or kill it (0, or just hit enter): '); if keepit==2 % Kill all casts subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out),bins_depth,'r-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out),bins_depth,'r-','linewidth',2); i_out_final=i_out; else % Ask user about each cast if(keepit) subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out(outloop)),bins_depth,'b.-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out(outloop)),bins_depth,'g.-','linewidth',2); else subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out(outloop)),bins_depth,'r-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out(outloop)),bins_depth,'r-','linewidth',2); i_out_final=[i_out_final i_out(outloop)]; end; for outloop=2:size(i_out,2) subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out(outloop)),bins_depth,'m-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out(outloop)),bins_depth,'m-','linewidth',2); keepit=input('Keep the magenta cast (1), or kill it (0, or just hit enter): '); if(keepit) subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out(outloop)),bins_depth,'b.-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out(outloop)),bins_depth,'g.-','linewidth',2); else subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out(outloop)),bins_depth,'r-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out(outloop)),bins_depth,'r-','linewidth',2); i_out_final=[i_out_final i_out(outloop)]; end; end; end; else % The case of just 1 cast to be decided keepit=input('Keep the magenta cast (1), kill it (0, or just hit enter): '); if(keepit) subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out(outloop)),bins_depth,'b.-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out(outloop)),bins_depth,'g.-','linewidth',2); else subplot(121); plot(allcasts(7:6+size(bins_depth,2),i_out(outloop)),bins_depth,'r-','linewidth',2); subplot(122); plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_out(outloop)),bins_depth,'r-','linewidth',2); i_out_final=[i_out_final i_out(outloop)]; end; end; % Build the final cast array, plot bad points on original map i_in=find(ismember([1:size(allcasts,2)],i_out_final)==0); qccasts=allcasts(:,i_in); num_elim_qc=size(allcasts,2)-size(qccasts,2); disp(' '); disp(['QC has eliminated ',num2str(num_elim_qc),' casts. ',num2str(size(qccasts,2)),' casts remain.']); else % None flagged i_out_final=[]; qccasts=allcasts; disp(' '); disp(['QC has not eliminated any casts. ',num2str(size(qccasts,2)),' casts remain.']); end; else % QC==0 - but print out casts if plotyes==1 qccasts=allcasts; disp(' '); disp(['Skipping quality control... ',num2str(size(qccasts,2)),' casts remain.']); if plotyes > 0 if exist('f_qc') figure(f_qc); clf; else f_qc=figure; end; subplot(121); hold on; plot(allcasts(7:6+size(bins_depth,2),:),bins_depth,'b.-'); title('Temperature casts'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); ylabel('Depth (m)'); xlabel('Temperature (C)'); subplot(122); hold on; plot(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),:),bins_depth,'b.-'); title('Salinity casts'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); ylabel('Depth (m)'); xlabel('Salinity'); % Compute mean profile and standard deviation and overlay on profile plots temp_mean=nanmean(allcasts(7:6+size(bins_depth,2),:),2); temp_std=nanstd(allcasts(7:6+size(bins_depth,2),:),0,2); figure(f_qc); subplot(121); plot(temp_mean,bins_depth,'k','linewidth',3); plot(temp_mean-temp_std,bins_depth,'k','linewidth',2); plot(temp_mean+temp_std,bins_depth,'k','linewidth',2); sal_mean=nanmean(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),:),2); sal_std=nanstd(allcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),:),0,2); figure(f_qc); subplot(122); plot(sal_mean,bins_depth,'k','linewidth',2); plot(sal_mean-sal_std,bins_depth,'k','linewidth',2); plot(sal_mean+sal_std,bins_depth,'k','linewidth',2); end; % plotyes end; % qc if plotyes==2 & qc print('-dpng',[saveprefix 'qc.png']); end; % Decimate num_casts=[]; for binloop=1:size(bins_dist,2) num_casts=[num_casts size(find(qccasts(1,:)==binloop),2)]; end; if(plotyes>0 & decim) f_decim=figure; bar('v6',bins_dist,num_casts); hold on; plot([-onshore offshore],[mean(num_casts) mean(num_casts)],'r'); title('Number of casts per cross-shelf bin (average number in red)'); ylabel('Casts'); xlabel('Horizontal bin distance (km)'); end; if(decim) % Decimation disp(' '); disp('You have chosen to decimate the data. The goal in this step is to reduce casts in horizontal bins with too many casts.'); disp(' '); disp('First, the mean number of casts for each bin will be computed. Then, the bins with more casts than this number will be sorted on julian date, and casts eliminated evenly in time to achieve the mean number of casts.'); disp(' '); disp(['The average number of casts per cross-shelf bin is ',num2str(mean(num_casts)),'.']); i_toomany=find(num_casts > mean(num_casts)); i_justright=find(num_casts <= mean(num_casts)); decimcasts=qccasts(:,find(ismember(qccasts(1,:),i_justright))); for binloop=1:size(i_toomany,2) % Extract all casts for each bin that needs to be decimated qccasts_temp=qccasts(:,find(ismember(qccasts(1,:),i_toomany(binloop)))); qccasts_temp_sorted=sortrows(qccasts_temp',6)'; % Find out how many to get rid of... keepeveryx=size(qccasts_temp_sorted,2) / floor(mean(num_casts)); indicestokeep=round(1:keepeveryx:size(qccasts_temp_sorted,2)); decimcasts=[decimcasts qccasts_temp_sorted(:,indicestokeep)]; end; num_casts_decim=[]; for binloop=1:size(bins_dist,2) num_casts_decim=[num_casts_decim size(find(decimcasts(1,:)==binloop),2)]; end; num_elim_decim=size(qccasts,2)-size(decimcasts,2); disp(' '); disp(['Decimation has eliminated ',num2str(num_elim_decim),' casts. ',num2str(size(decimcasts,2)),' casts remain.']); if plotyes > 0 figure(f_decim); bar('v6',bins_dist,num_casts_decim,.6,'y'); hold on; %plot([-onshore offshore],[mean(num_casts) mean(num_casts)],'r'); title('Decimated number of casts per cross-shelf bin (average number in red)'); ylabel('Casts'); xlabel('Horizontal bin distance (km)'); end; else decimcasts=qccasts; disp(' '); disp(['Skipping decimation... ',num2str(size(decimcasts,2)),' casts remain.']); end; % Decimation if plotyes==2 & decim print('-dpng',[saveprefix 'decimated.png']); end; % Plot, on the original map, in blue, the final casts that are going to be % used and the casts that were eliminated from QC in red if plotyes | length(bathychoice)==1 figure(f_map); m_plot(decimcasts(2,:),decimcasts(3,:),'b.'); if(qc & ~isempty(i_out_final)) m_plot(allcasts(2,i_out_final),allcasts(3,i_out_final),'r.'); end; title(['Months: ',season,'; points: BLUE=included, RED=QC-rejected, GRN=skipped']); end; if plotyes==2 print('-dpng',[saveprefix 'map.png']); end; % Now, finally, we have an array of CTD casts that is QCed and decimated % Compute the final mean and standard deviation for each bin disp(' '); disp('Computing means and standard deviations...'); [X,Y]=meshgrid(bins_dist,bins_depth); MEANTEMP=NaN*ones(size(bins_depth,2),size(bins_dist,2)); MEANSAL=MEANTEMP; MEANST=MEANTEMP; MEANC=MEANTEMP; STDTEMP=MEANTEMP; STDSAL=MEANTEMP; STDST=MEANTEMP; STDC=MEANTEMP; NUMPOINTS=MEANTEMP; for binloop=1:size(bins_dist,2) i_bin=find(decimcasts(1,:)==binloop); tempdata=decimcasts(7:6+size(bins_depth,2),i_bin); saldata=decimcasts(7+size(bins_depth,2):6+2*size(bins_depth,2),i_bin); densdata=decimcasts(7+2*size(bins_depth,2):6+3*size(bins_depth,2),i_bin); sspddata=decimcasts(7+3*size(bins_depth,2):6+4*size(bins_depth,2),i_bin); MEANTEMP(:,binloop)=nanmean(tempdata,2); STDTEMP(:,binloop)=nanstd(tempdata,0,2); MEANSAL(:,binloop)=nanmean(saldata,2); STDSAL(:,binloop)=nanstd(saldata,0,2); MEANST(:,binloop)=nanmean(densdata,2); STDST(:,binloop)=nanstd(densdata,0,2); MEANC(:,binloop)=nanmean(sspddata,2); STDC(:,binloop)=nanstd(sspddata,0,2); NUMPOINTS(:,binloop)=sum(~isnan(tempdata),2); end; % Assign NaNs to areas where NUMPOINTS < minnumpts MEANTEMP(NUMPOINTS < minnumpts)=NaN; STDTEMP(NUMPOINTS < minnumpts)=NaN; MEANSAL(NUMPOINTS < minnumpts)=NaN; STDSAL(NUMPOINTS < minnumpts)=NaN; MEANST(NUMPOINTS < minnumpts)=NaN; STDST(NUMPOINTS < minnumpts)=NaN; MEANC(NUMPOINTS < minnumpts)=NaN; STDC(NUMPOINTS < minnumpts)=NaN; blankfile=[AVGSEC(:,1:2); AVGSEC(length(AVGSEC),1) min(AVGSEC(:,2))-10; AVGSEC(1,1) min(AVGSEC(:,2))-10; AVGSEC(1,1:2)]; if sum(sum(isnan(MEANTEMP))) ~= size(bins_dist,2) * size(bins_depth,2) if smoothing > 0 save blankfile.temp blankfile -ascii eval(sprintf('GRID_ROPT=''-R%d/%d/%d/%d'';',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1))); eval(sprintf('GRID_IOPT=''-I%d/%d'';',hbinsize,vbinsize)); GRID_MOPT='-Mblankfile.temp'; eval(sprintf('SMOO_SOPT=''-S%d'';',smoothing)); mymat_t=[]; mymat_s=[]; mymat_d=[]; mymat_c=[]; mymat_tstd=[]; mymat_sstd=[]; mymat_dstd=[]; mymat_cstd=[]; for cloop=1:length(bins_dist) for rloop=1:length(bins_depth) if NUMPOINTS(rloop,cloop) >= minnumpts & isnan(MEANTEMP(rloop,cloop))==0 mymat_t=[mymat_t; X(rloop,cloop) Y(rloop,cloop) MEANTEMP(rloop,cloop)]; mymat_tstd=[mymat_tstd; X(rloop,cloop) Y(rloop,cloop) STDTEMP(rloop,cloop)]; mymat_s=[mymat_s; X(rloop,cloop) Y(rloop,cloop) MEANSAL(rloop,cloop)]; mymat_sstd=[mymat_sstd; X(rloop,cloop) Y(rloop,cloop) STDSAL(rloop,cloop)]; mymat_d=[mymat_d; X(rloop,cloop) Y(rloop,cloop) MEANST(rloop,cloop)]; mymat_dstd=[mymat_dstd; X(rloop,cloop) Y(rloop,cloop) STDST(rloop,cloop)]; mymat_c=[mymat_c; X(rloop,cloop) Y(rloop,cloop) MEANC(rloop,cloop)]; mymat_cstd=[mymat_cstd; X(rloop,cloop) Y(rloop,cloop) STDC(rloop,cloop)]; end; end; end; disp('Gridding and smoothing mean temperature using Matlab ppzgrid'); eval(sprintf('[Z,grd]=ppzinit(mymat_t,''-R%d/%d/%d/%d -I%d/%d -Mblankfile.temp'');',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1),hbinsize,vbinsize)); eval(sprintf('[ppz,grdout]=ppsmooth(Z,grd,''-S%d'');',smoothing)); MEANTEMP=flipud(ppz); disp('Gridding and smoothing std temperature using Matlab ppzgrid'); eval(sprintf('[Z,grd]=ppzinit(mymat_tstd,''-R%d/%d/%d/%d -I%d/%d -Mblankfile.temp'');',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1),hbinsize,vbinsize)); eval(sprintf('[ppz,grdout]=ppsmooth(Z,grd,''-S%d'');',smoothing)); STDTEMP=flipud(ppz); disp('Gridding and smoothing mean salinity using Matlab ppzgrid'); eval(sprintf('[Z,grd]=ppzinit(mymat_s,''-R%d/%d/%d/%d -I%d/%d -Mblankfile.temp'');',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1),hbinsize,vbinsize)); eval(sprintf('[ppz,grdout]=ppsmooth(Z,grd,''-S%d'');',smoothing)); MEANSAL=flipud(ppz); disp('Gridding and smoothing std salinity using Matlab ppzgrid'); eval(sprintf('[Z,grd]=ppzinit(mymat_sstd,''-R%d/%d/%d/%d -I%d/%d -Mblankfile.temp'');',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1),hbinsize,vbinsize)); eval(sprintf('[ppz,grdout]=ppsmooth(Z,grd,''-S%d'');',smoothing)); STDSAL=flipud(ppz); disp('Gridding and smoothing mean density using Matlab ppzgrid'); eval(sprintf('[Z,grd]=ppzinit(mymat_d,''-R%d/%d/%d/%d -I%d/%d -Mblankfile.temp'');',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1),hbinsize,vbinsize)); eval(sprintf('[ppz,grdout]=ppsmooth(Z,grd,''-S%d'');',smoothing)); MEANST=flipud(ppz); disp('Gridding and smoothing std density using Matlab ppzgrid'); eval(sprintf('[Z,grd]=ppzinit(mymat_dstd,''-R%d/%d/%d/%d -I%d/%d -Mblankfile.temp'');',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1),hbinsize,vbinsize)); eval(sprintf('[ppz,grdout]=ppsmooth(Z,grd,''-S%d'');',smoothing)); STDST=flipud(ppz); disp('Gridding and smoothing mean soundspeed using Matlab ppzgrid'); eval(sprintf('[Z,grd]=ppzinit(mymat_c,''-R%d/%d/%d/%d -I%d/%d -Mblankfile.temp'');',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1),hbinsize,vbinsize)); eval(sprintf('[ppz,grdout]=ppsmooth(Z,grd,''-S%d'');',smoothing)); MEANC=flipud(ppz); disp('Gridding and smoothing std soundspeed using Matlab ppzgrid'); eval(sprintf('[Z,grd]=ppzinit(mymat_cstd,''-R%d/%d/%d/%d -I%d/%d -Mblankfile.temp'');',bins_dist(1),bins_dist(end),bins_depth(end),bins_depth(1),hbinsize,vbinsize)); eval(sprintf('[ppz,grdout]=ppsmooth(Z,grd,''-S%d'');',smoothing)); STDC=flipud(ppz); end; if(plotyes>0) f1=figure; % Number of points per output grid node pcolor(X,Y,NUMPOINTS); shading flat; colorbar('eastoutside'); hold on; plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); ylabel('Depth (m)'); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Number of points for each bin']); if plotyes==2 print('-dpng',[saveprefix 'numpoints.png']); end; f2=figure; % Mean temperature [c,h]=contourf(X,Y,MEANTEMP,[floor(min(min(MEANTEMP))):1:ceil(max(max(MEANTEMP)))]); clabel(c,h); hold on; fill(blankfile(:,1),blankfile(:,2),[.7 .7 .7]) plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); ylabel('Depth (m)'); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Mean temperature (C)']); if plotyes==2 print('-dpng',[saveprefix 'meantemp.png']); end; f3=figure; % Standard deviation temperature [c,h]=contourf(X,Y,STDTEMP,[0:0.25:ceil(max(max(STDTEMP)))]); clabel(c,h); hold on; fill(blankfile(:,1),blankfile(:,2),[.7 .7 .7]) plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); ylabel('Depth (m)'); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Std dev temperature (C)']); if plotyes==2 print('-dpng',[saveprefix 'stdtemp.png']); end; f4=figure; % Mean salinity [c,h]=contourf(X,Y,MEANSAL,[floor(min(min(MEANSAL))):0.1:ceil(max(max(MEANSAL)))]); clabel(c,h); hold on; fill(blankfile(:,1),blankfile(:,2),[.7 .7 .7]) plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); ylabel('Depth (m)'); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Mean salinity']); if plotyes==2 print('-dpng',[saveprefix 'meansal.png']); end; f5=figure; % Standard deviation salinity [c,h]=contourf(X,Y,STDSAL,[0:0.05:ceil(max(max(STDSAL)))]); clabel(c,h); hold on; fill(blankfile(:,1),blankfile(:,2),[.7 .7 .7]) plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); ylabel('Depth (m)'); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Std dev salinity']); if plotyes==2 print('-dpng',[saveprefix 'stdsal.png']); end; f6=figure; % Mean sigma-t [c,h]=contourf(X,Y,MEANST-1000,[floor(min(min(MEANST-1000))):0.5:ceil(max(max(MEANST-1000)))]); clabel(c,h); hold on; fill(blankfile(:,1),blankfile(:,2),[.7 .7 .7]) plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); ylabel('Depth (m)'); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Mean sigma-t (kg/m^3)']); if plotyes==2 print('-dpng',[saveprefix 'meanst.png']); end; f7=figure; % Standard deviation sigma-t [c,h]=contourf(X,Y,STDST,[0:0.1:ceil(max(max(STDST)))]); clabel(c,h); hold on; fill(blankfile(:,1),blankfile(:,2),[.7 .7 .7]) plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); ylabel('Depth (m)'); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Std dev sigma-t (kg/m^3)']); if plotyes==2 print('-dpng',[saveprefix 'stdst.png']); end; f8=figure; % Mean sound speed [c,h]=contourf(X,Y,MEANC,[floor(min(min(MEANC))):5:ceil(max(max(MEANC)))]); clabel(c,h); hold on; fill(blankfile(:,1),blankfile(:,2),[.7 .7 .7]) plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); ylabel('Depth (m)'); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Mean sound speed (m/s)']); if plotyes==2 print('-dpng',[saveprefix 'meanc.png']); end; f9=figure; % Standard deviation sound speed [c,h]=contourf(X,Y,STDC,[0:1:ceil(max(max(STDC)))]); clabel(c,h); hold on; fill(blankfile(:,1),blankfile(:,2),[.7 .7 .7]) plot(AVGSEC(:,1),AVGSEC(:,2),'m'); ylab=get(gca,'ytick'); set(gca,'ytick',ylab); set(gca,'yticklabel',abs(ylab)); axis([-onshore offshore -abs(maxdepth) 0]); xlabel(['Distance from ',num2str(bathyline),'m isobath (km)']); ylabel('Depth (m)'); text(X(1,:),[ones(size(X(1,:)))*.75],'+'); title([filename,'; Months: ',season,'; Std dev sound speed (m/s)']); if plotyes==2 print('-dpng',[saveprefix 'stdc.png']); end; end; % if(plotyes>0) disp(' '); disp('Program finished.'); else disp(' '); disp('Program aborted: there are not enough casts to make a meaningful plot!'); end; disp(' '); disp(['Input CTD cast summary for ',filename]); disp(['Total extracted: ',num2str(totalextracted)]); disp([' Duplicates removed: ',num2str(num_elim_dupes)]); if(qc) disp([' QC removed: ',num2str(num_elim_qc)]); end; if(decim) disp([' Decimation removed: ',num2str(num_elim_decim)]); end; disp(['Total remaining: ',num2str(size(decimcasts,2))]); disp(' ');