How to reproject the data in Matlab?

17 次查看（过去 30 天）
Sophia
Sophia 2016-11-28
评论： Kelly Kearney 2016-12-16
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 个评论
Sophia
Sophia 2016-11-29
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Kelly Kearney
Kelly Kearney 2016-12-5
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 个评论
显示 1更早的评论
Sophia
Sophia 2016-12-14
编辑：Sophia 2016-12-14
*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 2016-12-16
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.

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2016-11-30
[Lat1, Long1] = meshgrid(lat1, long1);
[Lat2, Long2] = meshgrid(lat2, long2);
estimated_Data1 = interp2(Lat1, Long1, Data1, Lat2, Long2);
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2016-11-30
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 2016-12-5
编辑：Sophia 2016-12-5
Error using interp2/makegriddedinterp (line 222)
Input grid is not a valid MESHGRID.
Error in interp2 (line 133)
F = makegriddedinterp(X, Y, V, method);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Coordinate Reference Systems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by