grs80 = almanac('earth','grs80','km');
domeRadius =  300;  % km
domeLat =  20;       % degrees
domeLon = -77;       % degrees
domeAlt = 0;         % km

[x,y,z] = sphere(20);
xLV = domeRadius * x;
yLV = domeRadius * y;
zLV = domeRadius * z;
zLV(zLV < 0) = 0;
figure('Renderer','opengl')
surf(xLV, yLV, zLV,'FaceColor','yellow','FaceAlpha',0.5)
axis equal

[xECEF, yECEF, zECEF] = lv2ecef(xLV, yLV, zLV, ...
    domeLat * pi/180, domeLon * pi/180, domeAlt, grs80);
surf(xECEF, yECEF, zECEF,'FaceColor','yellow','FaceAlpha',0.5)
axis equal

figure('Renderer','opengl')
ax = axesm('globe','Geoid',grs80,'Grid','on', ...
    'GLineWidth',1,'GLineStyle','-',...
    'Gcolor',[0.9 0.9 0.1],'Galtitude',100);
set(ax,'Position',[0 0 1 1]);
axis equal off
view(3)

load topo
geoshow(topo,topolegend,'DisplayType','texturemap')
demcmap(topo)
land = shaperead('landareas','UseGeoCoords',true);
plotm([land.Lat],[land.Lon],'Color','black')
rivers = shaperead('worldrivers','UseGeoCoords',true);
plotm([rivers.Lat],[rivers.Lon],'Color','blue')

surf(xECEF, yECEF, zECEF,'FaceColor','yellow','FaceAlpha',0.5)