# Visualizing Geoid Height with given longitude, latitude and geoid heights

14 views (last 30 days)
Bora K on 5 Jan 2020
Answered: Meg Noah on 6 Jan 2020
Hello everyone,
I need to visualize the geoid height for a given txt file which has 260281x3 size. First column is longitude second is latitude and the third one is geoid heights. Could anyone help me to visualize the x as a longitude y as a latitude and the colour bar as geoid heights ? Thanks a lot in advance.
geoid = 'wwgrid.txt';
% Create 2-D plot using |meshm|
geoshow(geoid, 'DisplayType', 'texturemap');
%scatter(1:10,rand(1,10),50,rand(1,10),'*');

Meg Noah on 6 Jan 2020
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
% 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)