Home > atmlab > datasets > GriddedDataset.m

GriddedDataset

PURPOSE ^

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

GriddedDataset.m

SOURCE CODE ^

0001 classdef GriddedDataset < SatDataset
0002     % Class for gridded datasets
0003     %
0004     % This class contains functionality for gridded datasets.  Gridded
0005     % datasets are a special case of SatDatasets.  The special property of
0006     % a gridded dataset is that any datapoint can be uniquely identified by
0007     % lat, lon, time, where lat, lon, time are three cartesian axes.
0008 
0009     methods
0010         
0011         %% constructor
0012         
0013         function self = GriddedDataset(varargin)
0014             self = self@SatDataset(varargin{:});
0015         end
0016         
0017         %% overloaded methods
0018         
0019         function n = granule_first_line(varargin)
0020             n = int32(1);
0021         end        
0022         
0023         %% new methods
0024         
0025         function S = read_from_grid(self, lat, lon, time, fields)
0026             % Get data based on lat, lon, time.
0027             %
0028             % Currently, time is assumed to be on a resolution of one hour!
0029             %
0030             % FORMAT
0031             %
0032             %   S = ds.read_from_grid(lat, lon, time, fields)
0033             %
0034             % IN
0035             %
0036             %   lat     vector of latitudes
0037             %   lon     vector of longitudes, size must match lat
0038             %   time    vector or matrix of time: size must match lat&lon,
0039             %           and must be sorted.  Each row is a datevec.
0040             %   fields  Cell array of strings, fields to read
0041             %
0042             % OUT
0043             %
0044             %   S       structure with one member for each 'fields', each
0045             %           member is a vector with sizes matching lat/lon/time
0046             
0047             rqre_datatype(lat, @isvector);
0048             rqre_datatype(lon, @isvector);
0049             assert(length(lat)==length(lon), ['atmlab:' mfilename ':invalid'], 'lat and lon must be the same length, were found to differ');
0050             assert(size(time, 1)==length(lat), ['atmlab:' mfilename ':invalid'], 'no. rows in time must match length of lat and lon, were found to differ');
0051             assert(isequal(time, sortrows(time)), ['atmlab:' mfilename ':invalid'], 'time must be sorted in increasing order');
0052             
0053             % use time_plus_30m to get the nearest hour, for example,
0054             % between 23:30 and 24:00 nearest hour is on the following day
0055             time_plus_30m = datevec(datenum(time) + 1800/86400);
0056             [~, uni_indices] = unique(time_plus_30m(:, 1:3), 'rows');
0057             uni_indices = [uni_indices; length(lat)+1]; % add length(lat)+1 because later I run slices to llat-1
0058             S = struct();
0059             for i = 1:length(uni_indices)-1;
0060                 gran = self.find_granule_covering_instant(time_plus_30m(uni_indices(i), :), '');
0061                 if isempty(gran)
0062                     error(['atmlab:' mfilename ':nodata'], ...
0063                         ['No gridded data found for time %s %s. ' ...
0064                         'Please add to %s.'], ...
0065                         self.name, num2str(time_plus_30m(uni_indices(i), :)), ...
0066                         self.basedir);
0067                 end
0068                 SS = self.read_granule(gran, '', fields);
0069                 ustime = SS.epoch + SS.time;
0070                 assert(length(unique(sign(diff(SS.lat))))==1, ['atmlab:' mfilename ':invalid'], 'Weird lats came from source');
0071                 assert(length(unique(sign(diff(SS.lon))))==1, ['atmlab:' mfilename ':invalid'], 'Weird lons came from source');    
0072                 llat = lat(uni_indices(i):uni_indices(i+1)-1);
0073                 llon = lon(uni_indices(i):uni_indices(i+1)-1);
0074                 ttime = time(uni_indices(i):uni_indices(i+1)-1, :);
0075                 ttime_us = date2unixsecs(ttime(:, 1), ttime(:, 2), ttime(:, 3), ttime(:, 4), ttime(:, 5), ttime(:, 6));
0076                 llat_i = round(interp1(SS.lat, 1:length(SS.lat), llat, 'linear', 'extrap'));
0077                 % next line is an ugly hack, not 100% sure why this happens
0078                 llat_i(llat_i>length(SS.lat))=length(SS.lat);
0079                 llon_i = round(interp1(shift_longitudes(SS.lon, -180, 180), 1:length(SS.lon), shift_longitudes(llon, -180, 180), 'linear', 'extrap'));
0080                 ttime_i = round(interp1(ustime, 1:length(ustime), ttime_us, 'linear', 'extrap'));
0081                 if any(ttime_i<=0)
0082                     logtext(atmlab('ERR'), [...
0083                         'Warning: rounding %d timestamps in the first 30 minutes of the month up instead of down. ' ...
0084                         'This is because the first hour of a month is contained in the /previous/ month-file, ' ...
0085                         'and I haven''t fixed the sorting-per-day for those instances.\n'], sum(ttime_i<=0));
0086                     ttime_i(ttime_i<=0) = 1;
0087                 end
0088                 if any(llat_i==0)
0089                     logtext(atmlab('ERR'), [...
0090                          'Warning: found some latitudes so close to the' ...
0091                          'south pole that they''re out of range']);
0092                     llat_i(llat_i==0) = 1;
0093                 end
0094                 mask = true(size(llat_i));
0095                 mask(ttime_i>length(SS.time)) = false;
0096                 clear ii;
0097                 ii(mask) = sub2ind(size(SS.(fields{1})), llon_i(mask), llat_i(mask), ttime_i(mask));
0098                 ii(~mask) = NaN;
0099                 for j = 1:length(fields)
0100                     field = fields{j};
0101                     content(mask) = SS.(field)(ii(mask));
0102                     content(~mask) = NaN;
0103                     S.(field)(uni_indices(i):uni_indices(i+1)-1) = content;
0104                     clear content;
0105                     %S.(field)(uni_indices(i):uni_indices(i+1)-1) = SS.(field)(ii);
0106                 end
0107             end
0108         end
0109     end
0110 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005