First, create a 2D map of the heights were the rows are latitude and the columns are longitude.
You could use the matlab unique() on the list of latitudes and on the list of longitudes to create 1D arrays to define the row and column values. Then use the meshgrid command to make a 2D map of latitude and a 2D map of longitude. If the data are nicely ordered you can just use the reshape command to make the 1D array of heights into a 2D array. If it is not, then loop through lat and lon of the 2D array (each pixel) and search the Nx2 array accordingly doing a nearest neighbor interpolation into the 2D array of heights.
Most likely the data are nicely ordered so that reshape will work to create the 2D array. You might need to transpose the array. If each record first loops through a series of latitudes for one longitude, then increments longitude and loops through the same set of latitudes, then the height data. Below is some script code that creates fake topo data to illustrate:
nlong = 50;
nlat = 25;
[ICOL,IROW] = meshgrid(1:nlong,1:nlat);
longitude = zeros(nlat,nlong);
latitude = zeros(nlat,nlong);
height = zeros(nlat,nlong);
for ilong = 1:nlong
for ilat = 1:nlat
longitude(ilat,ilong) = 0.5*ilong - 70;
latitude(ilat,ilong) = 0.5*ilat + 40;
height(ilat,ilong) = 2i*pi*rand(1,1);
end
end
% making better fake data to illustrate
height = abs(ifft2(ifftshift((hypot(IROW,ICOL)+1e-5).^(-2.1).*exp(height))));
% shows ordering of longitude outer cycle and latitude inner cycle
% this is very likely how your data are organized
X = reshape(longitude,nlat*nlong,1);
Y = reshape(latitude,nlat*nlong,1);
Z = reshape(height,nlat*nlong,1);
for ipt = 1:nlat*nlong
fprintf(1,'%f %f %f\n',X(ipt),Y(ipt),Z(ipt));
end
% if this is how your data are arranged, then to create the 2D map:
latitude1D = unique(Y);
longitude1D = unique(X);
height2D = reshape(Z,nlat,nlong);
latitude2D = reshape(Z,nlat,nlong);
longitude2D = reshape(Z,nlat,nlong);
% if the reverse is true, that latitude is the outer cycle, then you need
% to transpose the result
% height2D = height2D';
%% referene matrix
% Reference matrix for mapping geoid heights to lat/lon on globe.
R1 = makerefmat('RasterSize',size(height2D), ...
'Latlim', [min(latitude(:)) max(latitude(:))], ...
'Lonlim', [min(longitude(:)) max(longitude(:))] );
%% visualization 1: DTED map with contour lines
figure();
hold on;
mapshow(zeros(size(height2D)),R1,'CData',height2D,'DisplayType','surface');
cW = mapshow(height2D,R1,'DisplayType','contour','linecolor','black', ...
'showtext','on');
demcmap(height2D);
axis equal
axis tight
%% visualization 2: with coastline
coast = load('coast');
% Create 2-D plot using |meshm|
h2D = figure;
set(h2D,'Position',[20 75 700 600],'Toolbar','figure');
gridDegInc = 0.5;
ast2DGeoidPlot(R1,height2D,coast,gridDegInc)
%% another visualization using the aerospace toolbox
figure
gcm = usamap([min(latitude1D) max(latitude1D)], [min(longitude1D) max(longitude1D)]);
hold on;
geoshow(gcm,height2D,R1,'DisplayType','texturemap');
demcmap(height2D);
daspectm('m',1);
geoshow(latitude2D,longitude2D,height2D,'DisplayType','contour', ...
'LineColor','black');
%% another visualization
figure
gcm = usamap([min(latitude1D) max(latitude1D)], [min(longitude1D) max(longitude1D)]);
hold on;
geoshow(gcm,height2D,R1,'Zdata',height2D, 'CData',height2D,'DisplayType','surface')
demcmap(height2D)
daspectm('m',1)
view(3)
Best Answer