Visualizing Geoid Height with given longitude, latitude and geoid heights

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';
load(geoid);
% Create 2-D plot using |meshm|
geoshow(geoid, 'DisplayType', 'texturemap');
%scatter(1:10,rand(1,10),50,rand(1,10),'*');

Answer 1

[Meg Noah]

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)