Figure 4: pandas / tabular data - scatter plot - seismicity around Indonesia

[yvonnefroehlich]

This PR adds a JN for pandas - tabular data - DataFrame:

• Figure.plot - scatter plot (color- and size-coding)
• Colorbar and legend for size-coding
• Figure.histogram (not only maps)

Preview:

[gitnotebooks]

Found 1 changed notebook. Review the changes at https://app.gitnotebooks.com/GenericMappingTools/pygmt-paper-figures/pull/5

[seisman]

@yvonnefroehlich Could you please also include the image in the notebooks, which can make review easier.

[yvonnefroehlich]

@yvonnefroehlich Could you please also include the image in the notebooks, which can make review easier.

I have included the output images within the notebooks for all PRs.

During preparing the AGU24 preworkshop we faced an issue with the review of JNs using gitnotebooks, please see GenericMappingTools/agu24workshop#9 (comment). It looks like that images lead to large diffs and for some JNs these diffs become too large and reviewing the JNs failed, but it worked when not showing the output images (less usage of fig.show or apply “clear all cells” before committing).

[seisman]

My thoughts on this example:

1. A real-world example is better than using a sample data. I suggest requesting a catalog from USGS, similar to what they Paul and others did in Figure 6 of the GMT6 paper). We can use the requests library to download the catalog from USGS on the fly, or we can just store the catalog in the repository, then load it into a pandas.DataFrame object.
2. We can request a catalog of mainshock and its aftershocks, then it would be a good opportunity to use the Figure.meca method to add a beachball of the mainshock.
3. I'm OK with adding a histogram in this example, and thinking if we should have a more complex layout instead? Below is a simple example, with an irregular grid. It can be set up via Figure.subplot, but maybe it's a good time to introduce Figure.shift_origin?

import pygmt

fig = pygmt.Figure()
fig.coast(projection="M10c", region=[131, 152, 33, 51], frame=True, land="gray")
# Manually tagging a subplot
fig.text(text="(a)", position="TL", font="12p,1", fill="white", offset="0.2/-0.2") 

fig.shift_origin(xshift="w+2c", yshift="h/2+0.8")
fig.basemap(region=[3, 8, 0, 40], projection="X5c/5c", frame=True)
fig.text(text="(b)", position="TL", font="12p,1", fill="white", offset="0.2/-0.2")

fig.shift_origin(yshift="-h-1.2")
fig.basemap(region=[0, 700, 0, 40], projection="X5c/5c", frame=True)
fig.text(text="(c)", position="TL", font="12p,1", fill="white", offset="0.2/-0.2")

fig.show()

[yvonnefroehlich]

1. A real-world example is better than using a sample data. I suggest requesting a catalog from USGS, similar to what they Paul and others did in Figure 6 of the GMT6 paper). We can use the requests library to download the catalog from USGS on the fly, or we can just store the catalog in the repository, then load it into a pandas.DataFrame object.

I also prefare to have real data. We also used real data for the seismology example in the workshop last year (wow - this is already nearly one year ago).

2. We can request a catalog of mainshock and its aftershocks, then it would be a good opportunity to use the Figure.meca method to add a beachball of the mainshock.

Sounds good.

3. I'm OK with adding a histogram in this example, and thinking if we should have a more complex layout instead? Below is a simple example, with an irregular grid. It can be set up via Figure.subplot, but maybe it's a good time to introduce Figure.shift_origin?

I agree that we should mention Figure.shift_origin in one of the examples, as this method can be quite convenient to create more complicated figure or subplot layouts. So, for this example, we have (i) a map with a scatterplot of the epicenters, (ii) and (iii) histograms for the moment magnitude and the hypocentral depth. Or do we only want one histogram and something else for the third panel?

[seisman]

1. Maybe we can also use Figure.inset here, depending on whether it makes the code too complicated.

Or do we only want one histogram and something else for the third panel?

Looking at GMT's plotting modules (https://docs.generic-mapping-tools.org/6.6/modules.html#plotting) and gallery examples (https://docs.generic-mapping-tools.org/6.6/gallery.html), I wish we could use GMT's rose module to show a rose diagram, but I don't know what kind of seismicity data is suitable for such diagrams.

3. I'm also looking at the images in your repository and find two figures: https://github.com/yvonnefroehlich/gmt-pygmt-plotting/tree/main/013_general_maps. I really hope that we could include them as our Figure 4 (well suited!) and Figure 7 (although your example doesn't use geopandas).

[yvonnefroehlich]

Or do we only want one histogram and something else for the third panel?

Looking at GMT's plotting modules (https://docs.generic-mapping-tools.org/6.6/modules.html#plotting) and gallery examples (https://docs.generic-mapping-tools.org/6.6/gallery.html), I wish we could use GMT's rose module to show a rose diagram, but I don't know what kind of seismicity data is suitable for such diagrams.

Hm, good questions. Not sure if there is something suitable for our example. For teleseismic studies, the backazimuth coverage by earthquakes around a recording station. But so far we plot the area with the epicenters.

3. I'm also looking at the images in your repository and find two figures: https://github.com/yvonnefroehlich/gmt-pygmt-plotting/tree/main/013_general_maps. I really hope that we could include them as our Figure 4 (well suited!) and Figure 7 (although your example doesn't use geopandas).

Wow, somebody is looking at these maps. Thanks for the positive feedback 🙂!
Just added the codes as they are in the scripts of the related PRs.

Regarding the epicenter map with the depth section at the sides: I saw such maps in publications and was wondering how people made the depth sections at the bottom and right. Thus I tried to figure out how to get a plot with a Cartesian depth axis but with geographical longitude or latitude axis to match the Mercator projection of the map in the center. As far as I know, there is no projection provided by GMT to do this directly, but you have to create several frames. Besides this, one needs to know the height of the map to know the height of the depth section on the right. Of course you can adjust this manually. But I think it could be nice to have a script where you can adjust the study area and the dimensions of the depth sections are adjusted automatically. GMT provides the mapproject module to get the width and height of a plot and save them into a file. But currently this module is not wrapped in PyGMT. Thus I call GMT within Python, create the file with the two values, and read the file. It works, but it's not really elegant :D. I think there is also something in the PyGMT codes with extracts the current plot dimensions.

Regarding the choropleth map: Here I actually tried to make a choropleth map by explicitly avoiding the polygon geometry provided by GeoPandas but using the dcw parameter from GMT (the script is very slow).

[seisman]

GMT provides the mapproject module to get the width and height of a plot and save them into a file. But currently this module is not wrapped in PyGMT.

mapproject is on the way! xref: GenericMappingTools/pygmt#4209, GenericMappingTools/pygmt#1986, GenericMappingTools/pygmt#4204.

[seisman]

I'm still debating which version to choose. The old one includes one scatter plot and two histograms. It not only demonstrates pandas support but also shows that PyGMT can generate histograms. The new one features three scatter plots and illustrates how to use the mapproject wrapper to obtain height and map coordinates, but it might be a bit too technical?

I'm a bit inclined toward the old version.

[yvonnefroehlich]

I'm still debating which version to choose. The old one includes one scatter plot and two histograms. It not only demonstrates pandas support but also shows that PyGMT can generate histograms. The new one features three scatter plots and illustrates how to use the mapproject wrapper to obtain height and map coordinates, but it might be a bit too technical?

I'm a bit inclined toward the old version.

I agree the advantage of the old version is that the histogram is a non-map plot. The new version is definitely more technical, but this kind of plot seems to be widely used in seismology publications. And getting the figure height can be useful for other figures too. Maybe we can include the mapproject wrapper in the old version: we can extract the height of the main map and use it to determine and set the height of the histograms (see commit d8238d7).

[yvonnefroehlich]

1. Maybe we can also use Figure.inset here, depending on whether it makes the code too complicated.

I added an inset map in commit 6c1bfad. It's only 4 lines more (without comments). But the inset has to be added before plotting the epicenters on the main map to have the dimensions of the main map available for Figure.shift_origin.

[seisman]

Can we just simplify it to:

Suggested change

	fig.colorbar(frame=["xaf+lHypocentral depth", "y+lkm"])
	fig.colorbar(frame="af+lHypocentral depth (km)")

[yvonnefroehlich]

Hm. I placed the unit as y label to avoid that we have to think about having units in brackets or using a /🙃😅.

[seisman]

OK, we can keep it now. But I feel the frame parameter is still not Pythonic and we may need to implement the Frame class anyway.

[yvonnefroehlich]

Oh, yes, the frame parameter is maybe the most complicated parameter in GMT / PyGMT and provides a huge amount of quite different functionalities 🙂.

[yvonnefroehlich]

Hm, maybe we should have a look at the author guidelines of $G^3$ if there are any rules how to set units in figures.

[seisman]

Do we really need to set the FONT here?

[yvonnefroehlich]

I would remove it; please see my comment #5 (comment). And in addition, there are different themes in GMT which also affect the defaults. I think having different font sizes look OK.

[seisman]

I feel this PR is good to merge except for two minor comments.

As before, please clean up the Python version and the output PNG images.

[seisman]

Line 51: "Magnitude M" seems duplicate, "Magnitude" is enough.