How to reproject the data in Matlab?

13 views (last 30 days)
Sophia
Sophia on 28 Nov 2016
Commented: Kelly Kearney on 16 Dec 2016
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.
  2 Comments
Sophia
Sophia on 29 Nov 2016
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

Sign in to comment.

Answers (2)

Kelly Kearney
Kelly Kearney on 5 Dec 2016
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.
  3 Comments
Sophia
Sophia on 14 Dec 2016
Edited: Sophia on 14 Dec 2016
*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
Kelly Kearney on 16 Dec 2016
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.

Sign in to comment.


Walter Roberson
Walter Roberson on 30 Nov 2016
[Lat1, Long1] = meshgrid(lat1, long1);
[Lat2, Long2] = meshgrid(lat2, long2);
estimated_Data1 = interp2(Lat1, Long1, Data1, Lat2, Long2);
  3 Comments
Walter Roberson
Walter Roberson on 30 Nov 2016
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
Sophia on 5 Dec 2016
Edited: Sophia on 5 Dec 2016
Error using interp2/makegriddedinterp (line 222)
Input grid is not a valid MESHGRID.
Error in interp2 (line 133)
F = makegriddedinterp(X, Y, V, method);

Sign in to comment.

Categories

Find more on Read, Write, and Modify Image in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!