Improvements in pvl_getISDdata.m, pvl_readtmy3.m, and pvl_spa.m

Correctly calculate distances between lat/lon points in pvl_getISDdata.m.
Allow pvl_readtmy3.m to read TMY3 files from before January 2015 and after January 2015.
Improve the default delta T calculation within pvl_spa.m

Author: DanRiley

pvl_getISDdata.m

@@ -1,19 +1,21 @@
-function outfilen = pvl_getISDdata(lat,lon,yr,archive)
+function outfilen = pvl_getISDdata(latitude,longitude,year,archive)
 %
 % getISDdata: fetch data for one year of integrated surface data for location
 %
 % Syntax
 %   data = getISDdata(latitude,longitude,year,archive)
 %
 % Inputs:
-%   latitude  - latitude of location of interest
-%   longitude - longitude of location of interest
+%   latitude  - latitude of location of interest in decimal degrees
+%   longitude - longitude of location of interest in decimal degrees
 %   year      - desired year
 %   archive   - string specifying path to local ISD archive.  Current 
 %               directory is used if archive is not specified.
 %
 % Output:
-%   outfilen - name for the returned data file
+%   outfilen - name for the returned data file in the form "USAF-WBAN-yyyy"
+%      where USAF is the USAF ID of the station and WBAN is the WBAN ID of 
+%      the system, as provided by the ISD.
 %
 % Notes:
 %   This function copies up to two files from NOAA's ftp site to the
@@ -40,18 +42,30 @@
 mod = vers(pos+6);
 if versyr<2013 || ((versyr==2013) && mod=='a')
     display(['Matlab ' vers ' detected'])
-    display('Function pvl_getISDdata requires R2013b or later')
-    return
+    error('Function pvl_getISDdata requires R2013b or later')
 end
-    
-%% Get index data
-% Change pwd to a local directory if you want to use an archive
-if nargin < 4
+%% Check the input data
+p = inputParser;
+p.addRequired('latitude',@(x) all(isscalar(x) & isnumeric(x) & x<=90 & x>=-90));
+p.addRequired('longitude', @(x) all(isscalar(x) & isnumeric(x) & x<=180 & x>=-180));
+p.addRequired('year',@(x) all(isscalar(x) & isnumeric(x) & (x >= 1901)));
+p.addOptional('archive', @(x) ischar(x));
+p.parse(latitude, longitude, year, varargin{:});
+
+defaultchecker = {'archive'};
+if ~any(strcmp(defaultchecker,p.UsingDefaults))
     archive = pwd;
 end
 
+lat = p.Results.latitude;
+lon = p.Results.longitude;
+yr = p.Results.year;
+
+%% Get index data
+% Change pwd to a local directory if you want to use an archive
+
 isdfile = 'isd-history.csv';
-fname = [archive '\' isdfile];
+fname = [archive filesep isdfile];
 if ~exist(fname,'file')
     % Get index file
     ftphandle = ftp('ftp.ncdc.noaa.gov');
@@ -69,17 +83,23 @@
 index.Properties.VariableNames{'ELEV_M_'} = 'ELEV';
 index.USAF = char(index.USAF);
 index.WBAN = char(index.WBAN);
-index.LAT(strcmp('',index.LAT)) = {'99999'};
-index.LON(strcmp('',index.LON)) = {'99999'};
-index.ELEV(strcmp('',index.ELEV)) = {'99999'};
+index.LAT(strcmp('',index.LAT)) = {'NaN'};
+index.LON(strcmp('',index.LON)) = {'NaN'};
+index.ELEV(strcmp('',index.ELEV)) = {'NaN'};
 index.LAT = str2num(char(index.LAT)); %#ok<*ST2NM>
 index.LON = str2num(char(index.LON));
 index.ELEV = str2num(char(index.ELEV));
 index.BEGIN = datenum(char(index.BEGIN),'yyyymmdd');
 index.END = datenum(char(index.END),'yyyymmdd');
 
-%% Find row for closest location (might iterate through several locations)
-dist = hypot(index.LAT-lat,index.LON-lon);
+%% Find row for closest location (might iterate through several locations) using haversine formula
+% calculate the argument for the haversine formula, make sure the argument
+% is between -1 and 1
+arg = sqrt((sind((index.LAT - lat) ./ 2)).^2 + cosd(index.LAT).*cosd(lat).*(sind((index.LON-lon)./2)).^2);
+arg(arg < -1) = -1;
+arg(arg > 1) = 1;
+% compute the distance in meters between the provided lat/lon and the each lat/lon from the isd.
+dist = 2 .* 6371008.8 .* asin(arg); % mean earth radius is 6371.0088 km
 [~,idx] = sort(dist);
 
 % Find first row that contains year (might iterate through several years)
@@ -99,10 +119,15 @@
 
 % 
 % get file from ftp server if not in archive
-fname = [archive '\' outfilen];
+fname = [archive filesep outfilen];
 if ~exist(fname,'file')
     addr = sprintf('ftp://ftp.ncdc.noaa.gov/pub/data/noaa/%4d/',yr);
-    gunzip([addr '/' outfilen '.gz'],archive);
+    try
+        gunzip([addr '/' outfilen '.gz'],archive);
+    catch
+        error(['Unable to reach the ftp site for the integrated surface database.' ...
+            ' The data file cannot be accessed.']);
+    end
     %% Alt approach if gunzip doesn't seem to be working
     % % cd(ftphandle,num2str(yr));
     % % gunzip(mget(ftphandle,[fname '.gz']));

pvl_readtmy3.m

@@ -15,6 +15,9 @@
 %   If a FileName is not provided, the user will be prompted to browse to
 %   an appropriate TMY3 file.
 %
+%   TMY3 format was changed slightly on January 19, 2015 ([2]) to include
+%   Present weather variables along with changes to other variables. 
+%
 %   Input
 %     FileName - an optional argument which allows the user to select which
 %     TMY3 format file should be read. A file path may also be necessary if
@@ -106,12 +109,23 @@
 %       TMYData.Lprecipquantity - The period of accumulation for the liquid precipitation depth field, hour
 %       TMYData.LprecipSource - See [1], Table 1-5, 8760x1 cell array of strings
 %       TMYData.LprecipUncertainty - See [1], Table 1-6
+%       TMYData.PresWth - Present weather code, see [2]. Only available in
+%         TMY3 files created after the January 2015 TMY3 update.
+%       TMYData.PresWthSource - Present weather code source, see [2]. Only
+%         available in TMY3 files created after the Jan. 2015 TMY3 update.
+%       TMYData.PresWthUncertainty - Present weather code uncertainty, see
+%         [2]. Only available in TMY3 files after the Jan. 2015 update.
 %
 % Reference
 % [1] Wilcox, S and Marion, W., 2008. Users Manual for TMY3 Data Sets, 
 % NREL/TP-581-43156, National Renewable Energy Laboratory.  Available at
 % <http://www.nrel.gov/docs/fy08osti/43156.pdf>.
 %
+% [2] National Solar Radiation Data Base, 1991-2005 Update: Typical
+% Meteorological Year 3. 
+% http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/tmy3/ retrieved June
+% 6, 2016.
+%
 %
 % See also
 %   DATEVEC  PVL_MAKELOCATIONSTRUCT  PVL_MAKETIMESTRUCT  PVL_READTMY2
@@ -136,15 +150,118 @@
 Header2 = textscan(FileID, '%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s',1,'Delimiter',',');
 DataLines = textscan(FileID, '%s%s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%f%s%f',8760,'Delimiter',',');
 ST = fclose(FileID);
-%% Read in the data from the struct created from the textscan
+
+% Parse out the header information since it is common to both formats
 TMYData.SiteID = Header1{1};
 TMYData.StationName = Header1{2};
 TMYData.StationState = Header1{3};
 TMYData.SiteTimeZone = Header1{4};
 TMYData.SiteLatitude = Header1{5};
 TMYData.SiteLongitude = Header1{6};
 TMYData.SiteElevation = Header1{7};
+%% Determine if the file is from pre-January 19, 2015 or post-20150119
+if isempty(strfind(char(DataLines{1,1}{1,1}), 'PresWth'))
+    % The string 'PresWth' was not found in DataLines{1,1}{1,1}, so we have a
+    % pre-2015 TMY3 file
+    TMYData = ParsePre2015TMY(DataLines);
+else
+    % The string 'PresWth' was found in DataLines{1,1}{1,1}, so we have a
+    % post-2015 TMY3 file. We'll need to re-import the file to accomodate
+    % the additional columns.
+    FileID = fopen(FilePathAndNameAndExt);
+    Header1 = textscan(FileID, '%f%q%q%f%f%f%f', 1,'Delimiter',',');
+    Header2 = textscan(FileID, '%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s',1,'Delimiter',',');
+    DataLines = textscan(FileID, '%s%s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%f%s%f%f%s%f',8760,'Delimiter',',');
+    ST = fclose(FileID);
+    
+    TMYData = ParsePost2015TMY(DataLines);
+end
+end
 
+function TMYData = ParsePre2015TMY(DataLines)
+
+%% Read in the data from the struct created from the textscan for TMY3 
+% This fuction will read pre-2015 TMY3 files and sort them into column
+% vectors. A different format is used for TMY3 files after the January 2015
+% update, and thus requires a different function.
+
+TMYData.DateString = DataLines{1};
+TMYData.TimeString = DataLines{2};
+TMYData.ETR = DataLines{3};
+TMYData.ETRN = DataLines{4};
+TMYData.GHI = DataLines{5};
+TMYData.GHISource = DataLines{6};
+TMYData.GHIUncertainty = DataLines{7};
+TMYData.DNI = DataLines{8};
+TMYData.DNISource = DataLines{9};
+TMYData.DNIUncertainty = DataLines{10};
+TMYData.DHI = DataLines{11};
+TMYData.DHISource = DataLines{12};
+TMYData.DHIUncertainty = DataLines{13};
+TMYData.GHillum = DataLines{14};
+TMYData.GHillumSource = DataLines{15};
+TMYData.GHillumUncertainty = DataLines{16};
+TMYData.DNillum = DataLines{17};
+TMYData.DNillumSource = DataLines{18};
+TMYData.DNillumUncertainty = DataLines{19};
+TMYData.DHillum = DataLines{20};
+TMYData.DHillumSource = DataLines{21};
+TMYData.DHillumUncertainty = DataLines{22};
+TMYData.Zenithlum = DataLines{23};
+TMYData.ZenithlumSource = DataLines{24};
+TMYData.ZenithlumUncertainty = DataLines{25};
+TMYData.TotCld = DataLines{26};
+TMYData.TotCldSource = DataLines{27};
+TMYData.TotCldUnertainty = DataLines{28};
+TMYData.OpqCld = DataLines{29};
+TMYData.OpqCldSource = DataLines{30};
+TMYData.OpqCldUncertainty = DataLines{31};
+TMYData.DryBulb = DataLines{32};
+TMYData.DryBulbSource = DataLines{33};
+TMYData.DryBulbUncertainty = DataLines{34};
+TMYData.DewPoint = DataLines{35};
+TMYData.DewPointSource = DataLines{36};
+TMYData.DewPointUncertainty = DataLines{37};
+TMYData.RHum = DataLines{38};
+TMYData.RHumSource = DataLines{39};
+TMYData.RHumUncertainty = DataLines{40};
+TMYData.Pressure = DataLines{41};
+TMYData.PressureSource = DataLines{42};
+TMYData.PressureUncertainty = DataLines{43};
+TMYData.Wdir = DataLines{44};
+TMYData.WdirSource = DataLines{45};
+TMYData.WdirUncertainty = DataLines{46};
+TMYData.Wspd = DataLines{47};
+TMYData.WspdSource = DataLines{48};
+TMYData.WspdUncertainty = DataLines{49};
+TMYData.Hvis = DataLines{50};
+TMYData.HvisSource = DataLines{51};
+TMYData.HvisUncertainty = DataLines{52};
+TMYData.CeilHgt = DataLines{53};
+TMYData.CeilHgtSource = DataLines{54};
+TMYData.CeilHgtUncertainty = DataLines{55};
+TMYData.Pwat = DataLines{56};
+TMYData.PwatSource = DataLines{57};
+TMYData.PwatUncertainty = DataLines{58};
+TMYData.AOD = DataLines{59};
+TMYData.AODSource = DataLines{60};
+TMYData.AODUncertainty = DataLines{61};
+TMYData.Alb = DataLines{62};
+TMYData.AlbSource = DataLines{63};
+TMYData.AlbUncertainty = DataLines{64};
+TMYData.Lprecipdepth = DataLines{65};
+TMYData.Lprecipquantity = DataLines{66};
+TMYData.LprecipSource = DataLines{67};
+TMYData.LprecipUncertainty = DataLines{68};
+%% Create a MATLAB datenum  from the string date and time
+TMYData.DateNumber = datenum(strcat(TMYData.DateString,'-',TMYData.TimeString),'mm/dd/yyyy-HH:MM'); %MATLAB datenum format
+end
+
+
+
+function TMYData = ParsePost2015TMY(DataLines)
+
+%% Read in the data from the struct created from the textscan
 TMYData.DateString = DataLines{1};
 TMYData.TimeString = DataLines{2};
 TMYData.ETR = DataLines{3};
@@ -213,5 +330,9 @@
 TMYData.Lprecipquantity = DataLines{66};
 TMYData.LprecipSource = DataLines{67};
 TMYData.LprecipUncertainty = DataLines{68};
+TMYData.PresWth = DataLines{69};
+TMYData.PresWthSource = DataLines{70};
+TMYData.PresWthUncertainty = DataLines{71};
 %% Create a MATLAB datenum  from the string date and time
-TMYData.DateNumber = datenum(strcat(TMYData.DateString,'-',TMYData.TimeString),'mm/dd/yyyy-HH:MM'); %MATLAB datenum format
+TMYData.DateNumber = datenum(strcat(TMYData.DateString,'-',TMYData.TimeString),'mm/dd/yyyy-HH:MM'); %MATLAB datenum format
+end