Plotting data above the Earth

I am trying to produce a figure similar to the one below, which I found in a publication:

I want to have a spherical Earth and then plot some data slightly above the surface. I have the data in both latitude,longitude as well as Cartesian X,Y,Z formats. I don’t really know where to start, as I could not find an example GMT script online which is similar to what I want. Any ideas are most welcome.

The data file is below.

output.txt (194.2 KB)

Here is my initial attempt:

#!/usr/bin/env bash


gmt begin ex52

  gmt pscoast ${proj} -B0 -Rd -W0.3p

  gmt psxy output.txt ${proj} -Rd -Gred -Sc0.1c -Wfaint

gmt end show

with the result:

Its unclear to me how to get the red line to go off the surface of Earth

#!/usr/bin/env bash


gmt begin ex52

  gmt pscoast ${proj} -B0 -Rd -W0.3p

  gmt psxy output.txt -JG80W/0N/12c -Rd -Gred -Sc0.1c -Wfaint -X-1c -Y-1c

gmt end show

Plotting the lines in a slightly larger map may work. Note that in the 2nd command, the map width is 12c and -X-1c/-Y-1c are also used to make sure that the centers of the two globes match.

The output is

Oooo, this is a good question! Seisman’s answer is kind of makes it look okay, but it’s not really plotting the data exactly where it would fall. This approach is like magnifying the data, relative to the map below it.

Maybe some mix of using your Orthographic projection in combination with feeding output.txt, in Cartesian coordinates, into psxyz could do it. Question to GMT authors is can psxyz handle -JG type of projections with lat,lon,height coordinates? I fear that psxyz won’t automatically mask data points that fall behind the visible Earth, because psxyz is totally unaware that the Earth is there.

It would seem that a pre-processing step might be required to determine whether a point (above the earth) is visible or not, relative to the view point. Not an insurmountable problem, but kind of an ugly one.

Excellent question.

I’ve not needed to make such a plot in a long time, but I’d be interested to learn how it can be done. Long ago, there were occasions when it was desired to display satellite orbits around a planet. For me, this was so long ago (pre-GMT) that I used MacDraw Pro to hand-make the diagram, approximately representing satellite altitudes and inclinations as best as I could. E.g., this is from a late 1980’s LAGEOS-II brochure (the size of the satellite is way-magnified):

Thank you for your answer. Unfortunately, the red line should wrap around the backside of Earth which it doesn’t in this example, but the data has all lat/lons in it

Just a quick thought: Maybe plotting that data on a custom larger ellipsoid might yield what you want? See PROJ_ELLIPSOID

Now possible from Julia.

using GMT

mat = gmtread("output.txt", i="2,3,4");
orbits(, lon0=-80, show=1)

1 Like

interesting thank you