Home > atmlab > deprecated > collocate.m

collocate

PURPOSE ^

collocate Return collocations between matrices

SYNOPSIS ^

function collocations = collocate(t1, lat1, long1, t2, lat2, long2, maxdist, maxtime)

DESCRIPTION ^

 collocate Return collocations between matrices

 collocate searches for collocations between the measurements in (t1, lat1,
 long1) and (t2, lat2, long2). Latitudes and longitudes should be in degrees,
 the time can be in any unit. 'collocate' considers footprints as points and
 defines a collocation according to a maximum distance ('maxdist', in km).
 The maximum time to consider a collocation is in 'maxtime' and should be in
 the same unit as t1 and t2.

 FORMAT collocations = collocate(t1, lat1, long1, t2, lat2, long2, ...
               maxdist, maxtime)

 OUT M       Nx4 matrix. Each row corresponds to a collocation, and has the
             rows and columns of the 1st and 2nd dataset respectively. 

 IN  t1      Array length L1 giving the times of the scanlines.
 IN  lat1    Matrix L1xW1 with latitude for L1 scanlines with W1 scans per
             line. The latitude should be in degrees.
 IN  long1   Matrix L1xW1 with longitude (-180, 180), L1 scans, W1 scans/line.
 IN  t2      Array of length L2 giving the times of the scanlines.
 IN  lat2    Matrix L2xW2 with latitude, L2 scans, W2 scans/line.
 IN  long2   Matrix L2xW2 with longitude (-180, 180), L2 scans, W2 scans/line.
 IN  maxdist Max distance (km) to consider a collocation.
 IN  maxtime Max time (units as for t1, t2) to consider a collocation.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

collocate.m

SOURCE CODE ^

0001 function collocations = collocate(t1, lat1, long1, t2, lat2, long2, maxdist, maxtime)
0002 
0003 % collocate Return collocations between matrices
0004 %
0005 % collocate searches for collocations between the measurements in (t1, lat1,
0006 % long1) and (t2, lat2, long2). Latitudes and longitudes should be in degrees,
0007 % the time can be in any unit. 'collocate' considers footprints as points and
0008 % defines a collocation according to a maximum distance ('maxdist', in km).
0009 % The maximum time to consider a collocation is in 'maxtime' and should be in
0010 % the same unit as t1 and t2.
0011 %
0012 % FORMAT collocations = collocate(t1, lat1, long1, t2, lat2, long2, ...
0013 %               maxdist, maxtime)
0014 %
0015 % OUT M       Nx4 matrix. Each row corresponds to a collocation, and has the
0016 %             rows and columns of the 1st and 2nd dataset respectively.
0017 %
0018 % IN  t1      Array length L1 giving the times of the scanlines.
0019 % IN  lat1    Matrix L1xW1 with latitude for L1 scanlines with W1 scans per
0020 %             line. The latitude should be in degrees.
0021 % IN  long1   Matrix L1xW1 with longitude (-180, 180), L1 scans, W1 scans/line.
0022 % IN  t2      Array of length L2 giving the times of the scanlines.
0023 % IN  lat2    Matrix L2xW2 with latitude, L2 scans, W2 scans/line.
0024 % IN  long2   Matrix L2xW2 with longitude (-180, 180), L2 scans, W2 scans/line.
0025 % IN  maxdist Max distance (km) to consider a collocation.
0026 % IN  maxtime Max time (units as for t1, t2) to consider a collocation.
0027 
0028 % Originally created by Gerrit Holl.
0029 % $Id$
0030 
0031 warning(['atmlab:' mfilename], 'old style function, being phased out, use OO way');
0032 
0033 
0034 %% data checks
0035 
0036 nrows1 = size(lat1, 1);
0037 ncols1 = size(lat1, 2);
0038 nrows2 = size(lat2, 1);
0039 ncols2 = size(lat2, 2);
0040 nel1 = numel(lat1);
0041 nel2 = numel(lat2);
0042 
0043 errid = 'atmlab:collocate';
0044 % check dimensions
0045 assert(size(long1, 1)==nrows1, errid, 'length long1 differs from length lat1');
0046 assert(size(long2, 1)==nrows2, errid, 'length long2 differs from length long1');
0047 assert(length(t1)==nrows1, errid, 'length t1 differs from length lat1');
0048 assert(length(t2)==nrows2, errid, 'length t2 differs from length long2');
0049 assert(size(long1, 2)==ncols1, errid, 'width long1 differs from width lat1');
0050 assert(size(long2, 2)==ncols2, errid, 'width long2 differs from width lat2');
0051 
0052 % check correct numbers
0053 assert(max(abs(long1(:)))<=180, errid, 'Invalid data in long1');
0054 assert(max(abs(lat1(:)))<=90, errid, 'Invalid data in lat1');
0055 assert(max(abs(long2(:)))<=180, errid, 'Invalid data in long2');
0056 assert(max(abs(lat2(:)))<=180, errid, 'Invalid data in lat2');
0057 
0058 %% make t, rows, cols at the same grid --> no need for ind2sub
0059 
0060 t1f = reshape(repmat(t1, [1 ncols1]), [1 nel1]);
0061 rows1f = reshape(repmat((1:nrows1)', [1 ncols1]), [1 nel1]);
0062 cols1f = reshape(repmat(1:ncols1, [nrows1 1]), [1 nel1]);
0063 
0064 t2f = reshape(repmat(t2, [1 ncols2]), [1 nel2]);
0065 rows2f = reshape(repmat((1:nrows2)', [1 ncols2]), [1 nel2]);
0066 cols2f = reshape(repmat(1:ncols2, [nrows2 1]), [1 nel2]);
0067 
0068 %% determine earth radius
0069 
0070 earth_radius = constants('EARTH_RADIUS')/1e3; % km
0071 
0072 %% find collocations
0073 
0074 % arbitrary start for pre-alloc, but getting more if needed later
0075 collocations = zeros(numel(lat1), 4, 'uint32');
0076 
0077 % Bin the measurements into bins, where the 'data' is the index of the
0078 % measurement. We only need to consider those bins where:
0079 % - instrument 1 has any measurements in the cell
0080 % - instrument 2 has any mearurements in the cell or a nearby cell
0081 % "Nearby" cell means a neighbouring cell, except near the poles, where it
0082 % can be much further away (because the cells are equirectangular)
0083 
0084 gridsize = 2;
0085 
0086 [grid1, lat_grid1, lon_grid1] = binning_fast(...
0087     struct(...
0088         'lat', lat1(:), ...
0089         'lon', long1(:), ...
0090         'data', uint32(1:numel(lat1)).', ...
0091         'gridsize', gridsize));
0092 n_grid_lats = length(lat_grid1);
0093 n_grid_lons = length(lon_grid1);
0094 grid2 = binning_fast(...
0095     struct(...
0096         'lat', lat2(:), ...
0097         'lon', long2(:), ...
0098         'data', uint32(1:numel(lat2)).', ...
0099         'gridsize', gridsize));
0100 
0101 % Check each cell where there is data for grid1 AND at least one
0102 % neighbouring cell, or the cell itself, has data for grid2
0103 
0104 count = 0;
0105 
0106 % the width and height of the cells as a function of latitude
0107 cell_width = 2 * pi * earth_radius * cosd(lat_grid1) / 360;
0108 cell_height = 2 * pi * earth_radius / 360;
0109 
0110 % how many cells of latitude to check for collocations
0111 cells_in_lat_range = ceil(maxdist ./ cell_width);
0112 c_lon = ceil(maxdist ./ cell_height);
0113 
0114 for i = 1:size(grid1, 1) % latitude
0115     c_lat = cells_in_lat_range(i);
0116     for j = 1:size(grid1, 2) % longitude
0117         if any(grid1{i, j})
0118             
0119             % find indices of points in this grid-cell
0120             in_grid1 = grid1{i, j};
0121             
0122             % which range of grid-cells to look for collocations?
0123             
0124             % for longitude (cols), depends on cells_in_range
0125             cols = j-c_lat:j+c_lat;
0126             cols = mod(cols, n_grid_lons);
0127             cols(cols==0) = n_grid_lons;
0128 
0129             % for latitude (rows), no dependence on lat
0130 
0131             rows = i-c_lon:i+c_lon;
0132             % edges do not wrap: if rows < 1
0133             rows(rows<1) = [];
0134             rows(rows>n_grid_lats) = [];
0135             
0136             % find indices of points in nearby grid-cells to consider
0137             in_grid2 = grid2(rows, cols);
0138             % flatten and convert to an array
0139             in_grid2 = vertcat(in_grid2{:});
0140             if isempty(in_grid2)
0141                 continue
0142             end           
0143             
0144             % find combinations of points that are near each other
0145             for p1 = in_grid1'
0146                 index1 = rows1f(p1);
0147                 colno1 = cols1f(p1);
0148                 shorttime = in_grid2(abs(t2f(in_grid2) - t1f(p1)) < maxtime);
0149                 near = shorttime(sphdist(lat1(p1), long1(p1), ...
0150                     lat2(shorttime), long2(shorttime), ...
0151                     earth_radius) < maxdist);
0152 
0153                 nnear = length(near);
0154                 if nnear
0155                     index2 = rows2f(near)';
0156                     colno2 = cols2f(near)';
0157                     collocations(count+1:count+nnear, :) = ...
0158                         [repmat(index1, [nnear 1]) repmat(colno1, [nnear 1]) ...
0159                         index2 colno2];
0160                     count = count + nnear;
0161                 end
0162             end % for each in grid1
0163         end % if we have any candidate-collocations in this cell
0164         if count > .8*size(collocations, 1) % pre-allocate more
0165             collocations = [collocations; ...
0166                 zeros(size(collocations, 1), 4, 'uint32')];
0167         end
0168     end % for all grid rows
0169 end % for all grid columns
0170 
0171 % removing 0's
0172 collocations(collocations(:, 1)==0, :)=[];

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