How to reproject the data in Matlab?

I have two different datasets for the same region, i want to reproject one of them according to the lat/long from the other file.

Data1 is 361*361, associated lat/long file of the size 361*361

Data 2 is 119*177, lat1/long1file of the 119*177

How can i reproject one of them according to the lat/long from other file.

[per isakson]

See How can I convert mapping data between different projections ...

[Sophia]

Sorry, the question should be how can i put them on common grid.

[X1,Y1] = meshgrid(lat1,lat);
[X2,Y2] = meshgrid(long1,long);

or

[X1,Y1] = meshgrid(lat1,long1);
[X2,Y2] = meshgrid(lat,long);

to use grid data later to bring them on the same grid

Answer 1

[Kelly Kearney]

What sort of spatial distribution do your datasets have? Is it a grid, but not one sitting along straight parallels and meridians?

Unfortunately, none of Matlab's interpolation routines are ideally suited for gridded-but-not-meshgrid/ndgrid-style datasets. If you have this type of data (and a lot of geographic data fits this description), then I think your best bet it to try scatteredInterpolant. Beware of problems along the edges of your grid, though; scatteredInterpolant doesn't like it when its input data includes a lot of colinear points (and this sort of grid will include a lot of colinearity), and can sometimes return very out-of-range values, even with a nearest neighbor interpolant.

[Kelly Kearney]

Also, keep in mind when using scatteredInterpolant that it assumes your x- and y-coordinates are equal, which is not the case with lat/lon coordinates. I usually convert to a cartesian coordinate system first (see geodetic2ecef function).

[Sophia]

*This is what i have, but there is an issue near zero longitude values*

Attached is the figure that show lat/lat1 and long/long1 in one plot

And the figure that show the error in the result

% change NSIDC longitude from 180 to -180 to match OSISAF
long(long==180)=-180;
% long  = long  +360;
% long1 = long1 +360;

% OSI-SAF output
UU = NaN.*ones([119 177]);
VV = NaN.*ones([119 177]);

% NSERC input
U = w_u_mean;
V = w_v_mean;

% alterantive 0-360 degrees longitude
long360 = long;
long360(long360<=0) = 360 + long360(long360<=0);

for i=2:length(lat1(:,1))-1,

      disp(i);

      for j=2:length(lat1(1,:))-1,

          % define area to look for NSIDC grid points
          lats = [lat1(i-1,j);  lat1(i+1,j);  lat1(i,j-1);  lat1(i,j+1)];
          lons = [long1(i-1,j); long1(i+1,j); long1(i,j-1); long1(i,j+1)];

          latsmax = max(lats);
          latsmin = min(lats);
          lonsmax = max(lons);
          lonsmin = min(lons);

          % tricky when we go from 180 to -180. So change longtitude
          % values to go from 0-360 for those cases
          if sign(lonsmax)~=sign(lonsmin),

              lons(lons<=0) = 360+lons(lons<=0); %alternative from 0-360 degrees

              lonsmax = max(lons);
              lonsmin = min(lons);

              % find the [row,col] combinations with the grid points
              % within the area bounds
              [row,col] = find(lat>latsmin & lat<latsmax & long360>lonsmin & long360<lonsmax);

          else
              % the normal situation when signs are not chaning
              [row,col] = find(lat>latsmin & lat<latsmax & long>lonsmin & long<lonsmax);

          end

          if isempty(col), % no values found

              UU(i,j) = NaN;
              VV(i,j) = NaN;

          else

              arc = []; % arclength in degrees
              Un = [];  % all the velocity values within the area
              Vn = [];

              % k number of values found
              for k = 1:length(col),
                  arc = [arc; distance(lat1(i,j),long1(i,j),lat(row(k),col(k)),long(row(k),col(k)))];
                  Un = [Un; U(row(k),col(k))];
                  Vn = [Vn; V(row(k),col(k))];
              end

              % calculate distance from the arclength
              dis = distdim(arc,'deg','km','earth');

              % inverse distance weighting
              UU(i,j) = sum(Un.*(1./dis.^2))./sum(1./dis.^2);
              VV(i,j) = sum(Vn.*(1./dis.^2))./sum(1./dis.^2);

          end

      end 

end

% magnitude
Mag = sqrt(UU.^2 + VV.^2);

%OSISAF-NSIDC
t = [] ;
for row = 1:size(m_wr,1);
    for col = 1:size(m_wr,2) ;
        if m_wr(row,col) == 0 ;
           t (row,col)= 0;
        else
           t(row,col) = m_wr(row,col)-Mag(row,col);
        end
    end
end

[Kelly Kearney]

It's a bit difficult to determine what's going wrong in your code, since we don't have the data you're using. However, here's a short example of how I'd go about an interpolation given the two data grids you include in the .fig file:

% Extract your data grids

hfig = open('~/Downloads/lat_long_nsidc_osisaf.fig')
h1 = findall(hfig, 'color', 'r');
h2 = findall(hfig, 'color', 'b');
A.lat = cell2mat(get(h1, 'ydata'));
A.lon = cell2mat(get(h1, 'xdata'));
B.lat = cell2mat(get(h2, 'ydata'));
B.lon = cell2mat(get(h2, 'xdata'));
close(hfig);

% Downsample B, just to save time in this example

B.lat = B.lat(1:10:end,1:10:end);
B.lon = B.lon(1:10:end,1:10:end);

% Add some fake data on the A grid

[nr,nc] = size(A.lon);
A.val = peaks(max(nr,nc));
A.val = A.val(1:nr,1:nc);

% Set up a function to calculate distance in geographic space, formatted
% for pdist2

ell = referenceEllipsoid('earth');
distfun = @(ltln1, ltln2) distance(ltln1(1), ltln1(2), ltln2(:,1), ltln2(:,2), ell);

% Project A data onto B grid, nearest neighbor

nb = numel(B.lat);
idx = zeros(nb,1);
nchunk = 100;
for ii = 1:(ceil(nb/nchunk))
    fprintf('Iteration %d of %d\n', ii, ceil(nb/nchunk));
    tmp = (1:nchunk)' + nchunk*(ii-1);
    tmp = tmp(tmp <= nb);
    [~, idx(tmp)] = pdist2([A.lat(:) A.lon(:)], [B.lat(tmp) B.lon(tmp)], distfun, 'Smallest', 1);
end

B.val = reshape(A.val(idx), size(B.lat));

% Plot 

Cst = load('coast');
figure('color', 'w');
ax(1) = axes('position', [0 0.5 1 0.5]);
axesm('ortho', 'origin', [90 0 0], 'frame', 'on');
pcolorm(A.lat, A.lon, A.val);
plotm(Cst.lat, Cst.long, 'k');

ax(2) = axes('position', [0 0 1 0.5]);
axesm('ortho', 'origin', [90 0 0], 'frame', 'on');
pcolorm(B.lat, B.lon, B.val);
plotm(Cst.lat, Cst.long, 'k');
set(ax, 'visible', 'off');

If you'd rather do inverse distance weighting as opposed to nearest neighbor, you can change the pdist2 calculation to return more points, as well as the distance values, and go from there.

The nearest neighbor calculation is slow... I usually break things into chunks (as I've done in this example), but even so, it's slow and memory-intensive. It could probably be sped up through use of quadtrees or the like, but there aren't any pre-written Matlab functions to do that in geographic space, so I've just suffered through the slowness of the exhaustive search myself.

Answer 2

[Walter Roberson]

[Lat1, Long1] = meshgrid(lat1, long1);
[Lat2, Long2] = meshgrid(lat2, long2);
estimated_Data1 = interp2(Lat1, Long1, Data1, Lat2, Long2);

[Sophia]

[Lat2, Long2] = meshgrid(lat2, long2);

It took forever to run this Error message, not enough space.

[Walter Roberson]

If your lat1 and lat2 and long1 and long2 are already 2 dimensional rather than vectors then just use

   result = interp2(lat1, long1, data, lat2, long2)

[Sophia]

Error using interp2/makegriddedinterp (line 222)
Input grid is not a valid MESHGRID.

Error in interp2 (line 133)
       F = makegriddedinterp(X, Y, V, method);