Making some Mars maps with pygmt (extended)¶

This tutorial page covers the basics of creating some maps and 3D plot of Mars (yes! Mars). The idea here is to demonstrate that you can use a simple sequence of commands with PyGMT, a Python wrapper for the Generic Mapping Tools (GMT), and with some public data about the topography of Mars, create your own maps, as well as compare this topography with what we know of our own planet Earth.

first, some options¶

You can run this notebook using your local pygmt installation, or via Binder, or even Google Colaboratory. See comments for each option below.

A) A short note if you are using COLAB

The version of python in COLAB is different from what the newer GMT needs to install along with pygmt. So, one way around this problem is to reinstall GMT from scratch, along with other important packages. This is done with this block of commands below.

comment out the first line of the block (%%script echo skipping) if you want to use colab

%%script echo skipping

# because I like to enjoy my coffee in silence, it takes time.  
# (3 runs averaged 6 minutes to install everything ! keep drinking your coffee)
# comment the %%capture line if you want to see the colab VM working
%%capture
!sudo apt update 
!sudo apt upgrade -y
!sudo apt install -y build-essential cmake libcurl4-gnutls-dev libnetcdf-dev gdal-bin libgdal-dev libfftw3-dev libpcre3-dev liblapack-dev libblas-dev libglib2.0-dev ghostscript ghostscript-x graphicsmagick ffmpeg xdg-utils
# clone gmt from source
!git clone --depth 50 https://github.com/GenericMappingTools/gmt
# cmake everything
!cmake /content/gmt
# build and install
!cmake --build . --target install

# and last but not least
!pip install pygmt

# and if you don't believe in it
!gmt --version
!python --version

skipping

# Also, if you are in colab or trying from your jupyter, you will need the Mars Topography (MOLA) already in Netcdf
# a copy of the original file distributed from the Mars Climate Database,
# from the European Space Agency under ESTEC contract 11369/95/NL/JG(SC) and Centre National D'Etude Spatial
# is on GitHub.

!wget https://github.com/andrebelem/PlanetaryMaps/raw/v1.0/mola32.nc

B) Now, if you are using Binder or in your local jupyter

You just skip the block above. Make sure you have the mola32.nc in your folder.

%matplotlib inline

Now, how about Hawaii?¶

When we read about Olympus Mons, it is usually compared to Everest here on Earth. However, the most interesting thing is to compare it with another mountain range taking as a reference the abyssal seabed (without the ocean) - Hawaii. Interestingly, in terms of latitudes and longitudes on the planet, these two features are in almost the same position. To match the approximate dimensions, let’s crop a sample of the Earth Global Relief using pygmt.datasets with slices of 12 degrees of latitude and 30 degrees of longitude.

# get SRTM around Hawaii 
topo_hawaii = pygmt.datasets.load_earth_relief(region=[-170,-140,13,25],resolution="05m")

# and get the whole Earth at the same resolution of our low resolution Mars dataset
topo_globe = pygmt.datasets.load_earth_relief(region=[-180,180,-90,90],resolution="01d")

grdblend [NOTICE]: Remote data courtesy of GMT data server oceania [http://oceania.generic-mapping-tools.org]

grdblend [NOTICE]: SRTM15 Earth Relief v2.6 at 05x05 arc minutes reduced by Gaussian Cartesian filtering (26.2 km fullwidth) [Tozer et al., 2019].

grdblend [NOTICE]:   -> Download 180x180 degree grid tile (earth_relief_05m_p): S90W180

grdcut [NOTICE]: Remote data courtesy of GMT data server oceania [http://oceania.generic-mapping-tools.org]

grdcut [NOTICE]: SRTM15 Earth Relief v2.6 at 1x1 arc degrees reduced by Gaussian Cartesian filtering (314.5 km fullwidth) [Tozer et al., 2019].

grdcut [NOTICE]:   -> Download grid file [112K]: earth_relief_01d_p.grd

And we use the same sequence as above to make a map.

# second things, second

fig = pygmt.Figure()
# note I can add projection, after cmap and after, frame (and control frame)
fig.grdimage(grid=topo_hawaii,projection='M12c',frame='a5f1',cmap='geo')
# also, I can add a colorbar (later)
fig.colorbar(frame=["a2500", "x+lElevation", "y+lm"])

fig.show()

Another few notes

As we want to make a comparison, let’s keep the same color scale as Mars, still using as a basis for the Z plane, -5000 meters (see the line plane="-5000+ggrey" exactly like the map above. The inset in the upper right corner is the same and we adjust the bounding box coordinates used to cut out the topography, drawing in red on the map.

fig = pygmt.Figure()
# note I can add projection, after cmap and after, frame (and control frame)
frame =  ["xa5f1","ya5f1", "z5000+lmeters", "wSEnZ"]

topo_cpt = pygmt.makecpt(cmap='sealand',series=f'-5000/24000/1000',continuous=True)

fig.grdview(grid=topo_hawaii,
            region=[-170,-140,13,25,-5000,23000],
            frame=frame,
            perspective=[150,45],
            projection='M15c',
            zsize='4c',
            surftype='s',
            cmap=topo_cpt,
            plane="-5000+ggrey",
            shading='+a100+nt1',
            # Set the contour pen thickness to "0.1p"
            contourpen="0.1p",)

fig.colorbar(perspective=True, frame=["a5000", "x+lElevation", "y+lm"])

bounds = [[-170.,13.],
          [-170.,25.],
          [-140.,25.],
          [-140.,13.],
          [-170.,13.]]

with fig.inset(position="JTR+w3.5c+o0.2c", margin=0, box=None):
    # Create a figure in the inset using the global projection centered at Olympus MOns
    fig.grdimage(grid=topo_globe,region='g',frame='g',projection='G-160/19/3.5c"')
    fig.coast(region='g',shorelines="thin", frame="g")
    fig.plot(bounds,pen="1p,red")
    
fig.show()

Combining the two maps side by side¶

Basically it’s the same blocks as above, just using pygmt’s Figure.set_panel mechanism to tile.

fig = pygmt.Figure()

with fig.subplot(
    nrows=1, ncols=2, figsize=("28c", "16c"), autolabel=True, margins="1c"
):
    with fig.set_panel(panel=0):

        topo_cpt = pygmt.makecpt(cmap='sealand',series=f'-5000/24000/1000',continuous=True)

        frame =  ["xa5f1","ya5f1", "z5000+lmeters", "wSEnZ"]

        fig.grdview(grid=dset_olympus,
                    region=[210,240,13,25,-5000,23000],
                    frame=frame,
                    perspective=[150,45],
                    projection='M',
                    zsize='4c',
                    surftype='s',
                    cmap=topo_cpt,
                    plane="-5000+ggrey",
                    shading='+a100+nt1',
                    # Set the contour pen thickness to "0.1p"
                    contourpen="0.1p",)

        # we don't need the colormap in both figures
        #fig.colorbar(perspective=True, frame=["a5000", "x+lElevation", "y+lm"])
        
    with fig.set_panel(panel=1):
        frame =  ["xa5f1","ya5f1", "z5000+lmeters", "wSEnZ"]

        topo_cpt = pygmt.makecpt(cmap='sealand',series=f'-5000/24000/1000',continuous=True)

        fig.grdview(grid=topo_hawaii,
                    region=[-170,-140,13,25,-5000,23000],
                    frame=frame,
                    perspective=[150,45],
                    projection='M',
                    zsize='4c',
                    surftype='s',
                    cmap=topo_cpt,
                    plane="-5000+ggrey",
                    shading='+a100+nt1',
                    # Set the contour pen thickness to "0.1p"
                    contourpen="0.1p",)

        fig.colorbar(perspective=True, frame=["a5000", "x+lElevation", "y+lm"])       

fig.show()

Bonus map¶

Recently the rover Zhurong from the Tianwen-1’s mission landed successfully at 109.926°E, 25.066°N, in southern Utopia Planitia on Mars (check out the article of Ye, B., Qian, Y., Xiao, L., Michalski, J. R., Li, Y., Wu, B., & Qiao, L. (2021). Geomorphologic exploration targets at the Zhurong landing site in the southern Utopia Planitia of Mars. Earth and Planetary Science Letters, 576, 117199. https://doi.org/10.1016/j.epsl.2021.117199). We can create a map of the region with the landing point.

First, let’s locate Utopia Planitia. Take a look at Figure 1 by Ye et al. (2021).

fig = pygmt.Figure()

# we are using a Orthographic view centered at the landing site
fig.grdimage(grid=dset_mars_topo,region='g',frame='g',projection='G109.926/25.066/12c"',shading='+a100+nt1')

zhurong = [109.926,25.066]
Olympus = [360-210,19.0] #position for Olympus Mons - see the letf border of the area

# and we drop a "star" in the landing site and write with a small displacement of text
fig.plot(x=zhurong[0],y=zhurong[1],style="a0.5c", pen="1p,black", color="darkorange")
fig.text(x=zhurong[0]+5,y=zhurong[1]+5, text="Zhurong", font='10p,Helvetica-Bold')

fig.text(x=Olympus[0],y=Olympus[1], text="Olympus Mons", font='10p,Helvetica-Bold')

fig.colorbar(frame=["a5000", "x+lElevation", "y+lm"])
fig.show()

Crafting beautiful maps with PyGMT

Making some Mars maps with pygmt (extended)

Contents

Making some Mars maps with pygmt (extended)¶

first, some options¶

Mars dataset¶

Now, how about Hawaii?¶

Combining the two maps side by side¶

Bonus map¶

additional maps¶