Tutorial for maps of Mars using pygmt

book/_toc.yml

@@ -12,4 +12,7 @@ parts:
   chapters:
     - file: first-figure
     - file: ecosystem
+    - file: mars_maps
+      sections:
+        - file: mars_maps_extended
     - file: lidar_to_surface

book/mars_maps.ipynb

@@ -0,0 +1,207 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "19f4dff3",
+   "metadata": {},
+   "source": [
+    "# Making some Mars maps with pygmt\n",
+    "\n",
+    "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."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b48dd300",
+   "metadata": {},
+   "source": [
+    "## Mars dataset\n",
+    "\n",
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2f913d8d",
+   "metadata": {
+    "tags": [
+     "remove-stdout"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# Also, if you are in colab or trying from your jupyter, you will need the Mars Topography (MOLA) already in Netcdf\n",
+    "# a copy of the original file distributed from the Mars Climate Database,\n",
+    "# from the European Space Agency under ESTEC contract 11369/95/NL/JG(SC) and Centre National D'Etude Spatial\n",
+    "# is in the gdrive.\n",
+    "\n",
+    "!gdown 1fDzz8AxR1T58y0IGPhmbb1ZwrTLckp2G"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11b04ce1",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "import xarray as xr\n",
+    "\n",
+    "dset_mars = xr.open_dataset('mola32.nc')\n",
+    "dset_mars"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e9d6ef48",
+   "metadata": {},
+   "source": [
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8e7374cf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pygmt\n",
+    "\n",
+    "dset_mars_topo = dset_mars.alt.astype(float)\n",
+    "dset_mars_topo = pygmt.grdsample(grid=dset_mars_topo, translate=True, spacing=[1,1])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "700d1d1d",
+   "metadata": {},
+   "source": [
+    "Here we can create a map of the entire Martian surface, in the same projections we use for our planet."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "feecaa82",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "fig = pygmt.Figure()\n",
+    "\n",
+    "fig.grdimage(grid=dset_mars_topo,region='g',frame=True,projection='W12c')\n",
+    "fig.colorbar(frame=['a5000','x+lElevation','y+lm'])\n",
+    "\n",
+    "fig.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3e3530ac",
+   "metadata": {},
+   "source": [
+    "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`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "01fac891",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dset_olympus = dset_mars.sel(latitude=slice(25,13),longitude=slice(210,240)).alt.astype(float)\n",
+    "dset_olympus.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fe62c4a7",
+   "metadata": {},
+   "source": [
+    "We use the same sequence as other pygmt tutorial notebooks to make a map."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9f18b522",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = pygmt.Figure()\n",
+    "\n",
+    "fig.grdview(grid=dset_olympus,\n",
+    "           region=[210,240,13,25,-5000,23000],\n",
+    "           surftype='s',\n",
+    "           projection='M12c',\n",
+    "           perspective=[150,45],\n",
+    "           zsize='4c',\n",
+    "           frame=['xa5f1','ya5f1','z5000+lmeters','wSEnZ'],shading='+a50+nt1')\n",
+    "\n",
+    "bounds = [[210,13],\n",
+    "         [210,25],\n",
+    "         [240,25],\n",
+    "         [240,13],\n",
+    "         [210,13]]\n",
+    "\n",
+    "fig.colorbar(perspective=True,frame=['a5000','x+lElevation','y+lm'])\n",
+    "with fig.inset(position='JTR+w3.5c+o0.2c',margin=0,box=None):\n",
+    "    fig.grdimage(grid=dset_mars_topo,region='g',frame=True,projection='W3.5c')\n",
+    "    fig.plot(bounds, pen='1p,red')\n",
+    "\n",
+    "fig.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "55a14a22",
+   "metadata": {},
+   "source": [
+    "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.\n",
+    "\n",
+    "**A few notes**\n",
+    "\n",
+    "`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.\n",
+    "\n",
+    "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."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ac07a1c0",
+   "metadata": {},
+   "source": [
+    "There are a few bonus maps in the *extended version* in github.\n",
+    "\n",
+    "Have fun."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}