Given a matrix of LATITUTDE, LONGITUDE, and ELEVATION data for a planet, how can I make a 3d mesh out of it?

The above screenshot is a snippet of a matrix,"undo-index", whose size can be as small as 1x3 or as large as 5,000,000x3. Each row represents data points on Mar's surface. The first column is the latitude (values range from -90 to +90), second is longitude (values range from -180 to +180), and third is elevation/height from a predetermined "zero-line". My question is: what is the best way I could render all of these points in 3d-space and create a mesh out of them? I'm wanting it to look like this. I'm not sure how the creator of that picture made it but that's exactly what I want.

I'm also having the issue of using scatter3() to plot the points without the mesh surface. The below code is meant to convert the latitude,longitude angles and elevation into x,y,z cartesian coordinates. The problem is, the scatter3 plot is not right.

latitude = undo_index(:,1);
longitude = undo_index(:,2);
radius = undo_index(:,3);

x = cosd(latitude) .* cosd(longitude) .* radius;
y = cosd(latitude) .* sind(longitude) .* radius ;
z = sind(latitude) .* radius;

scatter3(x,y,z);

I'm not sure why the plot looks like that... it does not look like Mars terrain at all. I'm not sure if I have the equations incorrect or if I somehow assigned the wrong elevation to each latitude,longitude coordinate - which would explain why the plot looks so wrong. The thing is, I don't see how the problem could be caused by the latter because when I put the data points on a mercador projection, the map looks like the below, which is correct.

If I could get some help on how to solve my scatter3 problem as well as how I could add on the aforementioned mesh, that would be very much appreciated.

I've attached my commented code along with the original data set as a Google Drive Link from which I created the "undo_index" matrix - if anyone would like to take a look at that. Thank you in advance for any help and happy holidays!

FYI, calculate_coord.m is a Function file that I call from in main.m

Answer 1

[Mathieu NOE]

hello

I used this (excellent) Fex submission to plot your data : Surface Fitting using gridfit - File Exchange - MATLAB Central (mathworks.com)

what we want is to plot radius vs lon and lat so you can skip / delete these lines

% x = cosd(latitude) .* cosd(longitude) .* radius;
% y = cosd(latitude) .* sind(longitude) .* radius ;
% z = sind(latitude) .* radius;

Result obtained so far :

Code updated :

clc;
clear;
close all;
%Constants
elevations = load('elevation.mat').elevation;
dlat =1.013895028871900;
dlon =-0.081103695129408;
resolution = 10;%determines angle spacing between points in degreees. For example, points can be spaced 10 degrees from each other laterally and horizontally
grid_size = [floor(180/resolution), floor(360/resolution)];

%initials
lat = 0;
lon = -180;
coord_data = zeros(length(elevations), 3);

%Pairs elevations with the corresponding lat,lon coordinate in coord_data
for ii=1:length(elevations)
    coord_data(ii,:) = [lat,lon,elevations(ii)];
    [lat, lon, dlat, dlon] = calculate_coord(lat, lon, dlat, dlon);
end

%Rounds the lat,lon coordinates to the nearest multiple of resolution
coord_data(:,[1 2]) = resolution*round(coord_data(:,[1 2])/resolution);

%maps lat,lon coordinates to indices
ind_coord_data(:,1) = floor((90 + coord_data(:,1))/resolution)+1;
ind_coord_data(:,2) = floor((180 + coord_data(:,2))/resolution)+1;
ind_coord_data(:,3) = coord_data(:,3);

%Any invalid indices (aka row/col == 0) gets turned to 1
ind_coord_data(ind_coord_data(:,[1 2])<=0) = 1;

%If a index has multiple values, it gets assigned to the mean of
%all of its values
accum_coord_data = accumarray(ind_coord_data(:,[1 2]), ind_coord_data(:,3), [], @mean);

% imagesc(-180:180, -90:90,accum_coord_data)
% set(gca, 'YDir', 'normal', 'XTick', -180:30:180, 'YTick', -90:30:90);
% colormap(turbo);
% grid on;

%converts index back into lat,lon coordinates
[row, col, val] = find(accum_coord_data);
undo_index = [row, col, val];
undo_index(:,1) = resolution * (undo_index(:,1)-1) - 90;
undo_index(:,2) = resolution * (undo_index(:,2)-1) - 180;


%conversion of lat,lon,elevation to cartesian x,y,z coordinates for
%plotting

latitude = undo_index(:,1);
longitude = undo_index(:,2);
radius = undo_index(:,3);

% x = cosd(latitude) .* cosd(longitude) .* radius;
% y = cosd(latitude) .* sind(longitude) .* radius ;
% z = sind(latitude) .* radius;

% scatter3(x,y,z);
% axis equal

figure(1)
scatter3(longitude,latitude,radius,100,radius,'filled');
view(2)

% Turn the scatter point data into a surface
x = longitude;
y = latitude;
z = radius;
gx = linspace(min(x),max(x),100);
gy = linspace(min(y),max(y),100);
tic
g=gridfit(x,y,z,gx,gy);
figure(2)
surf(gx,gy,g);
axis([-180 180 -90 90]);
xlabel('Longitude');
ylabel('Latitude');
view(2)
camlight right;
lighting phong;
shading interp
line(x,y,z,'marker','.','markersize',4,'linestyle','none');
title 'MARS Surface'
hcb=colorbar('ver'); % colorbar handle
hcb.FontSize = 10;
hcb.Title.String = "Elevation (m)";
hcb.Title.FontSize = 12;

[Mathieu NOE]

I didn't really pay too much attention to the main code (the lat / lon data generation) but I was surprised to see that the raw data (elevations , size = 5000000x1) would at the end only result in 703 valid points . Why so much loss ?