vis/ts_at.m

   1 function ts = ts_at(lon, lat, xlon, xlat, v)
   2     % TS_AT interpolates time series of a variable from a grid to a point.
   3     %
   4     % This function performs interpolation of a variable defined on a 2D grid
   5     % to a specific geographic point (longitude, latitude) over a series of
   6     % time steps.
   7     %
   8     % Arguments:
   9     %   lon (1D array) - The longitude of the location to interpolate to.
  10     %   lat (double) - The latitude of the location to interpolate to.
  11     %   xlon (2D double array) - The longitudes of the grid points.
  12     %   xlat (2D double array) - The latitudes of the grid points.
  13     %   v (3D double array) - The values on the grid. The third dimension
  14     %                         is time.
  15     %
  16     % Returns:
  17     %   ts (1D double array) - The interpolated time series of the variable
  18     %                           at the specified location.
  19     %
  20     % Example:
  21     %   % Given a set of grid points (xlon, xlat) and variable values v(time)
  22     %   lon = -120.5; lat = 34.5;
  23     %   ts_values = ts_at(lon, lat, grid_lon, grid_lat, variable_values);
  24
  25     % Check input dimensions
  26     if ndims(v) ~= 3
  27         error('The variable v must be a 3D array where the last dimension is time.');
  28     end
  29
  30     % Find the four nearest grid points
  31     [numRows, numCols] = size(xlon);
  32     [~, idx] = min(abs(xlon(:) - lon) + abs(xlat(:) - lat));
  33     [iRow, iCol] = ind2sub([numRows, numCols], idx);
  34
  35     % Define the square subgrid indices
  36     iRow = max(iRow-1, 1):min(iRow+1, numRows);
  37     iCol = max(iCol-1, 1):min(iCol+1, numCols);
  38
  39     % Extract the subgrid coordinates
  40     subXlon = xlon(iRow, iCol);
  41     subXlat = xlat(iRow, iCol);
  42
  43     % Initialize the output time series
  44     numTimeSteps = size(v, 3);
  45     ts = zeros(numTimeSteps, 1);
  46
  47     % Loop over each time step
  48     for t = 1:numTimeSteps
  49         % Extract the values for the current time step
  50         subV = v(iRow, iCol, t);
  51
  52         % Flatten the subgrid values to column vectors
  53         subV = subV(:);
  54
  55         % Create an interpolant for the current time step
  56         F = scatteredInterpolant(subXlon(:), subXlat(:), subV);
  57
  58         % Perform the interpolation for the current time step
  59         ts(t) = F(lon, lat);
  60     end
  61 end