Making some Mars maps with pygmt
Contents
Making some Mars maps with pygmt¶
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.
Mars dataset¶
First, we open the mola32.nc file using xarray. Note the longitudes are from 0-360°, latitudes are distributed from North to South and the alt variable is the MOLA Topography at 32 pixels/degree built from original MOLA file megt90n000fb.img. The source is the Mars Climate Database from the European Space Agency.
# 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 in the gdrive.
!gdown 1fDzz8AxR1T58y0IGPhmbb1ZwrTLckp2G
import xarray as xr
dset_mars = xr.open_dataset('mola32.nc')
dset_mars
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/file_manager.py:210, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
209 try:
--> 210 file = self._cache[self._key]
211 except KeyError:
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/lru_cache.py:56, in LRUCache.__getitem__(self, key)
55 with self._lock:
---> 56 value = self._cache[key]
57 self._cache.move_to_end(key)
KeyError: [<class 'netCDF4._netCDF4.Dataset'>, ('/home/runner/work/egu22pygmt/egu22pygmt/book/mola32.nc',), 'r', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False)), 'a35790ca-f6de-41bd-af70-6130f4dd6705']
During handling of the above exception, another exception occurred:
FileNotFoundError Traceback (most recent call last)
Cell In[2], line 3
1 import xarray as xr
----> 3 dset_mars = xr.open_dataset('mola32.nc')
4 dset_mars
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/api.py:526, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, inline_array, backend_kwargs, **kwargs)
514 decoders = _resolve_decoders_kwargs(
515 decode_cf,
516 open_backend_dataset_parameters=backend.open_dataset_parameters,
(...)
522 decode_coords=decode_coords,
523 )
525 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 526 backend_ds = backend.open_dataset(
527 filename_or_obj,
528 drop_variables=drop_variables,
529 **decoders,
530 **kwargs,
531 )
532 ds = _dataset_from_backend_dataset(
533 backend_ds,
534 filename_or_obj,
(...)
542 **kwargs,
543 )
544 return ds
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/netCDF4_.py:577, in NetCDF4BackendEntrypoint.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, format, clobber, diskless, persist, lock, autoclose)
557 def open_dataset(
558 self,
559 filename_or_obj,
(...)
574 autoclose=False,
575 ):
576 filename_or_obj = _normalize_path(filename_or_obj)
--> 577 store = NetCDF4DataStore.open(
578 filename_or_obj,
579 mode=mode,
580 format=format,
581 group=group,
582 clobber=clobber,
583 diskless=diskless,
584 persist=persist,
585 lock=lock,
586 autoclose=autoclose,
587 )
589 store_entrypoint = StoreBackendEntrypoint()
590 with close_on_error(store):
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/netCDF4_.py:382, in NetCDF4DataStore.open(cls, filename, mode, format, group, clobber, diskless, persist, lock, lock_maker, autoclose)
376 kwargs = dict(
377 clobber=clobber, diskless=diskless, persist=persist, format=format
378 )
379 manager = CachingFileManager(
380 netCDF4.Dataset, filename, mode=mode, kwargs=kwargs
381 )
--> 382 return cls(manager, group=group, mode=mode, lock=lock, autoclose=autoclose)
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/netCDF4_.py:329, in NetCDF4DataStore.__init__(self, manager, group, mode, lock, autoclose)
327 self._group = group
328 self._mode = mode
--> 329 self.format = self.ds.data_model
330 self._filename = self.ds.filepath()
331 self.is_remote = is_remote_uri(self._filename)
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/netCDF4_.py:391, in NetCDF4DataStore.ds(self)
389 @property
390 def ds(self):
--> 391 return self._acquire()
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/netCDF4_.py:385, in NetCDF4DataStore._acquire(self, needs_lock)
384 def _acquire(self, needs_lock=True):
--> 385 with self._manager.acquire_context(needs_lock) as root:
386 ds = _nc4_require_group(root, self._group, self._mode)
387 return ds
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/contextlib.py:119, in _GeneratorContextManager.__enter__(self)
117 del self.args, self.kwds, self.func
118 try:
--> 119 return next(self.gen)
120 except StopIteration:
121 raise RuntimeError("generator didn't yield") from None
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/file_manager.py:198, in CachingFileManager.acquire_context(self, needs_lock)
195 @contextlib.contextmanager
196 def acquire_context(self, needs_lock=True):
197 """Context manager for acquiring a file."""
--> 198 file, cached = self._acquire_with_cache_info(needs_lock)
199 try:
200 yield file
File /usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/xarray/backends/file_manager.py:216, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
214 kwargs = kwargs.copy()
215 kwargs["mode"] = self._mode
--> 216 file = self._opener(*self._args, **kwargs)
217 if self._mode == "w":
218 # ensure file doesn't get overridden when opened again
219 self._mode = "a"
File src/netCDF4/_netCDF4.pyx:2353, in netCDF4._netCDF4.Dataset.__init__()
File src/netCDF4/_netCDF4.pyx:1963, in netCDF4._netCDF4._ensure_nc_success()
FileNotFoundError: [Errno 2] No such file or directory: b'/home/runner/work/egu22pygmt/egu22pygmt/book/mola32.nc'
Just like any other notebook with pygmt, we import the library and manipulate other data. To make a map of the entire Martian surface without a lot of time and memory, let’s reduce the resolution using grdsample. We also take the opportunity to transform an alt variable into a float to be used in maps.
import pygmt
dset_mars_topo = dset_mars.alt.astype(float)
dset_mars_topo = pygmt.grdsample(grid=dset_mars_topo, translate=True, spacing=[1,1])
Here we can create a map of the entire Martian surface, in the same projections we use for our planet.
fig = pygmt.Figure()
fig.grdimage(grid=dset_mars_topo,region='g',frame=True,projection='W12c')
fig.colorbar(frame=['a5000','x+lElevation','y+lm'])
fig.show()
A very interesting feature is Mount Olympus (Olympus Mons - see more details at https://mars.nasa.gov/resources/22587/olympus-mons), centered at approximately 19°N and 133°W, with a height of 22 km (14 miles) and approximately 700 km (435 miles) in diameter. Let’s use the original dataset at 32 pixels/degree resolution and plot a (not so interesting) map with xarray.
dset_olympus = dset_mars.sel(latitude=slice(25,13),longitude=slice(210,240)).alt.astype(float)
dset_olympus.plot()
We use the same sequence as other pygmt tutorial notebooks to make a map.
fig = pygmt.Figure()
fig.grdview(grid=dset_olympus,
region=[210,240,13,25,-5000,23000],
surftype='s',
projection='M12c',
perspective=[150,45],
zsize='4c',
frame=['xa5f1','ya5f1','z5000+lmeters','wSEnZ'],shading='+a50+nt1')
bounds = [[210,13],
[210,25],
[240,25],
[240,13],
[210,13]]
fig.colorbar(perspective=True,frame=['a5000','x+lElevation','y+lm'])
with fig.inset(position='JTR+w3.5c+o0.2c',margin=0,box=None):
fig.grdimage(grid=dset_mars_topo,region='g',frame=True,projection='W3.5c')
fig.plot(bounds, pen='1p,red')
fig.show()
And we’re going to add some perspective, as well as a more interesting color scale. For ease of understanding, let’s separate the region of interest with the same cutout that we created the base of the Olympus Mons topography dataset.
A few notes
zsize is a bit critical here because the volcano is very big (28 km if we consider -5000 to +23000 m). Likewise, perspective=[150,45] was chosen attempting (it’s a matter of taste) and depends of which flank of the volcano you want to show. But this choice has to be made according to shading since to give a good 3D impression, the lighting must be adjusted according to the elevation and azimuth of the perspective. Finally, the pen outline is made smooth and small to enhance the contours of the topography.
Finally, let’s make a combined map showing the planet in an inset in the upper right corner. We use the same bounding box coordinates used to cut out the topography, drawing in red on the map. Obviously here the color scale is the same.
There are a few bonus maps in the extended version in github.
Have fun.