Tutorial 4 - Geophysics (Seismology): Add JN

[yvonnefroehlich]

This PR adds the Jupyter notebook for Tutorial 4 - Geophysics (Seismology).

• Study region: Chilean subduction zone since there are significant topographic variation, earthquakes, and volcanoes
• Orientation: https://www.pygmt.org/dev/gallery/images/cross_section.html

Content:

• Prepare gridded data: pygmt.xyz2grd
• Create a contour map: pygmt.Figure.grdcontour
• Create a profile plot: pygmt.project and pygmt.grdtrack
• Additional features: seismicity, beachball, plate motion

Feel free to fill the subsections with suggestions for code and docs🙂.

Preview:

• https://github.com/GenericMappingTools/agu24workshop/blob/add-jn-geo-seis/book/tut4_geophysics.ipynb
• https://gitnotebooks.com/GenericMappingTools/agu24workshop/pull/8

Fixes: #6

[jhtong33]

Hi @yvonnefroehlich
do we use grdproject or project ?
Since grdproject is specified for a region rather than individual point-to-point.

see also: https://www.pygmt.org/v0.13.0/gallery/images/cross_section.html#sphx-glr-gallery-images-cross-section-py

[yvonnefroehlich]

Hi @yvonnefroehlich do we use grdproject or project ? Since grdproject is specified for a region rather than individual point-to-point.

see also: https://www.pygmt.org/v0.13.0/gallery/images/cross_section.html#sphx-glr-gallery-images-cross-section-py

Yep, I am sorry, I mean project I mistyped this 🙁.

[yvonnefroehlich]

Thanks a lot, @jhtong33 , for getting started here 👍!

I made a few updates:

• I prefer to use the Scientific colour maps by F. Crameri, i.e., for topography, "oleron".
• I reduced the resolution of the image showing within the Jupyter notebook to keep the file size smaller.
• I figured out that list input is supported for the center and endpoint parameters of project.
• Fix typo grdproject to project.
• So minor code style improvements.
• I think we should be consistent and use either single (') or double (") quotation marks; I prefer double quotation marks.

[jhtong33]

Thanks a lot, @jhtong33 , for getting started here 👍!

I made a few updates:

• I prefer to use the Scientific colour maps by F. Crameri, i.e., for topography, "oleron".
• I reduced the resolution of the image showing within the Jupyter notebook to keep the file size smaller.
• I figured out that list input is supported for the center and endpoint parameters of project.
• Fix typo grdproject to project.
• So minor code style improvements.
• I think we should be consistent and use either single (') or double (") quotation marks; I prefer double quotation marks.

Great! Thanks for your updates.
I used a remote dataset to convert data into a .txt file, which I then provided as input materials for pygmt.xyz2grd.

[jhtong33]

Hi @yvonnefroehlich
I have one idea that plotting beachballs for magnitude >6 events to make the tutorial more comprehensive.
However, If time make it challenging to cover all details within a 45-minute session, it would be great to focus on one significant event as an example.

see also: https://www.isc.ac.uk/cgi-bin/web-db-run?request=COMPREHENSIVE&out_format=CATCSV&searchshape=RECT&bot_lat=-40&top_lat=0&left_lon=-90&right_lon=-60&ctr_lat=&ctr_lon=&radius=&max_dist_units=deg&srn=&grn=&start_year=2022&start_month=1&start_day=01&start_time=00%3A00%3A00&end_year=2023&end_month=1&end_day=01&end_time=00%3A00%3A00&min_dep=&max_dep=&min_mag=6&max_mag=&req_mag_type=MW&req_mag_agcy=&include_links=on

[yvonnefroehlich]

Hm. But for me, this approach does not work in the expected way. Actually, this reminds me a bit of ObsPy where users actively need to copy the Stream object, before applying any processing steps in case they want to keep their raw data in the initial variable unchanged (https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.copy.html).

I misunderstood your meaning previously. I also encounter same condition in Jupyter lab. In obspy, they use copy.deepcopy(self). I tried to use same way to solve this issue, but the result is same.

This seem to be more complicated to implemente, please see GenericMappingTools/pygmt#3645 (comment) for details.

[yvonnefroehlich]

Are there any major review requests on this tutorial 🙂? Otherwise, I feel we can move on and merge this tutorial in the current state and, if needed, open separat PRs for minor corrections and improvements.

[Esteban82]

I have two observations only:

1. The first row of the seismicity_2022.csv file is: time,lon,lat,depth,mag. With GMT, it is good to add a # to the beginning of the line so that it is automatically identified as a header. I suspect that with PyGMT this is not necessary. is it true?
2. It would be possible to add a shading effect to the topographic map of South America. The idea would be to show that with GMT you can make very good maps simply.
   Then, during the workshop it might be better to say that one removes the shading effect to make the processing faster.

[yvonnefroehlich]

1. The first row of the seismicity_2022.csv file is: time,lon,lat,depth,mag. With GMT, it is good to add a # to the beginning of the line so that it is automatically identified as a header. I suspect that with PyGMT this is not necessary. is it true?

In this tutorial, we read the CSV file via pandas into a DataFrame. The first line is by default considered as the header line and used for the column names. Later, these column names can be used to access the single columns. pandas does not require a # to indicate a header line, but would include the # into the first column name 😬:

                            #time      lon      lat    depth  mag
0     2022-12-31T23:04:12.784000Z -69.6800 -16.5378  196.743  4.9
1     2022-12-31T12:47:57.842000Z -71.5625 -32.2686   48.300  4.5
2     2022-12-31T03:31:43.830000Z -68.8969 -23.1036   96.224  5.1
3     2022-12-31T00:43:57.999000Z -77.7617  -1.7451  182.543  4.3
4     2022-12-30T22:35:25.866000Z -67.1397 -24.0267  232.941  4.5
                          ...      ...      ...      ...  ...

Thus, I prefer to leave this file as is.

[yvonnefroehlich]

2. It would be possible to add a shading effect to the topographic map of South America. The idea would be to show that with GMT you can make very good maps simply.
   Then, during the workshop it might be better to say that one removes the shading effect to make the processing faster.

Hm. You mean something simliar as Max does in the Xarray tutorial in the section "Processing raster data" (https://github.com/GenericMappingTools/agu24workshop/blob/add-xarray-tutorial/book/tut03_spec_xarray.ipynb) using grdgradient? Maybe we can include this in the section "Additional comments" as the tutorial is already a bit longer?

[Esteban82]

Hm. You mean something simliar as Max does in the Xarray tutorial in the section "Processing raster data" (https://github.com/GenericMappingTools/agu24workshop/blob/add-xarray-tutorial/book/tut03_spec_xarray.ipynb) using grdgradient? Maybe we can include this in the section "Additional comments" as the tutorial is already a bit longer?

My idea is just to make the maps look better. I don't think it's necessary to explain anything. So, I think you can add something like shading="+a" in the code so the map looks better.

[yvonnefroehlich]

Hm. You mean something simliar as Max does in the Xarray tutorial in the section "Processing raster data" (https://github.com/GenericMappingTools/agu24workshop/blob/add-xarray-tutorial/book/tut03_spec_xarray.ipynb) using grdgradient? Maybe we can include this in the section "Additional comments" as the tutorial is already a bit longer?

My idea is just to make the maps look better. I don't think it's necessary to explain anything. So, I think you can add something like shading="+a" in the code so the map looks better.

Ah, I first thought about grdgradient but applying the illumination directly via grdimage sounds good to me. Made an update in commit a3ded9d.

[yvonnefroehlich]

I am wondering if and how it is possible to add links to the other Jupyter notebooks or scripts (if this is possible, maybe we can also do this in a common PR).

[jhtong33]

I'm not sure I fully understand your point. Are you suggesting that this information could be linked to the current gallery?

[yvonnefroehlich]

Hm. What do you mean with "current gallery"?

We partly mention other tutorials, which is good. I feel it would be nice to directly insert the link. But we have the URL from GitHub (https://github.com/GenericMappingTools/agu24workshop/blob/main/book/tut0XYZ.ipynb) as well as the one from the Jupyter book (https://www.generic-mapping-tools.org/agu24workshop/XYZ.html). So depending on which website people are, the link does lead them to another "base" website. As the Jupyter book website is the "official" website for this workshop and GitHub is more for development, maybe using the URL to the Jupyter book website makes more sense here?

[jhtong33]

Now I understand what you mean.

As the Jupyter book website is the "official" website for this workshop and GitHub is more for development, maybe using the URL to the Jupyter book website makes more sense here?

YES. Providing the URL to the Jupyter Book website will make it more convenient and readable.

[weiji14]

Wonderful, thanks @yvonnefroehlich and @jhtong33 for working on this! Sorry for leaving this last, Seismology isn't my area of expertise, so I was hoping someone would jump in to review, but I trust that this is ok. Personally, I would have preferred it if the .txt and .csv file were not in the repository since they are big files, but hopefully everyone has good internet connection during the workshop 😆

I've got some minor suggestions below, which either of you can apply, and then I will try to merge this by tomorrow. After this, we can work on polishing all of the notebooks to have a consistent header, title and such.

[jhtong33]

Personally, I would have preferred it if the .txt and .csv file were not in the repository since they are big files, but hopefully everyone has good internet connection during the workshop 😆

If needed, I can adjust the period of the earthquake catalog, which might reduce the file size. (current size: 63 kb)

[yvonnefroehlich]

Wonderful, thanks @yvonnefroehlich and @jhtong33 for working on this! Sorry for leaving this last, Seismology isn't my area of expertise, so I was hoping someone would jump in to review, but I trust that this is ok.

Now worries @weiji14, @Esteban82 made already some helpful suggestions! And somethings it's actually also good to have a review from another field and also good to a have native-speaking person reading through the documentation.

Personally, I would have preferred it if the .txt and .csv file were not in the repository since they are big files, but hopefully everyone has good internet connection during the workshop 😆

Yeah! Also was thinking about what is better: people download the data during the tutorials or we provide data files within the repo, when creating the pands / geopandas tutorial. Downloading the data directly underlines that PyGMT integrates smoothly in the complete analysis workflow, and we keep the repo size small. Regarding a stable internet connection, I feel (hope) this should work (I have never been to AGU or any other big conference yet, but a working (acceptable fast) internet connection is essential for a nicely running conference) .

Regarding the elevation data: The "normal" way doing this would be to simply download the built-in GMT remote dataset for the Earth relief into a xarray.DataArray and using it for the further processing and analysis. However, in this tutorial, we like to show how people can convert any specific grid which is provided as xyz columns into a [Py]GMT-ready grid. Thus, the elevation data was written to a txt file with xyz columns, and in the tutorial we convert it back to a [Py]GMT-ready grid via xyz2grd.

Regarding the seismicity data: ObsPy is required to request and download this data. If we include this part within the tutorial, ObsPy will become a dependency, which people need to install. This was probably the reason why Jing-Hui downloaded the data externally and wrote it to a CSV file for usage in the tutorial (reading via pandas).

I've got some minor suggestions below, which either of you can apply, and then I will try to merge this by tomorrow. After this, we can work on polishing all of the notebooks to have a consistent header, title and such.

Sounds good!

[jhtong33]

Regarding the elevation data: The "normal" way doing this would be to simply download the built-in GMT remote dataset for the Earth relief into a xarray.DataArray and using it for the further processing and analysis. However, in this tutorial, we like to show how people can convert any specific grid which is provided as xyz columns into a [Py]GMT-ready grid. Thus, the elevation data was written to a txt file with xyz columns, and in the tutorial we convert it back to a [Py]GMT-ready grid via xyz2grd.

Yvonne mentioned a critical point: many current datasets (e.g., crustal thickness and seismic tomography) are available only as xyz columns rather than pre-processed GMT grids. If beginners fail to correctly convert xyz data into a grid, using grdimage will not produce accurate results. That’s why providing the xyz file is essential.

Regarding the seismicity data: ObsPy is required to request and download this data. If we include this part within the tutorial, ObsPy will become a dependency, which people need to install. This was probably the reason why Jin-Hui downloaded the data externally and wrote it to a CSV file for usage in the tutorial (reading via pandas).

If including Obspy, installing and downloading IRIS dataset more relies on internet.
I can mitigate this by reducing the file size, perhaps by limiting the earthquake data to a shorter time period.

[weiji14]

Regarding the elevation data: The "normal" way doing this would be to simply download the built-in GMT remote dataset for the Earth relief into a xarray.DataArray and using it for the further processing and analysis. However, in this tutorial, we like to show how people can convert any specific grid which is provided as xyz columns into a [Py]GMT-ready grid. Thus, the elevation data was written to a txt file with xyz columns, and in the tutorial we convert it back to a [Py]GMT-ready grid via xyz2grd.

Yvonne mentioned a critical point: many current datasets (e.g., crustal thickness and seismic tomography) are available only as xyz columns rather than pre-processed GMT grids. If beginners fail to correctly convert xyz data into a grid, using grdimage will not produce accurate results. That’s why providing the xyz file is essential.

Ah I see, then teaching xyz2grd would make sense here. Thanks for the extra context, this is why having people working in the same field makes sense 😉

Regarding the seismicity data: ObsPy is required to request and download this data. If we include this part within the tutorial, ObsPy will become a dependency, which people need to install. This was probably the reason why Jin-Hui downloaded the data externally and wrote it to a CSV file for usage in the tutorial (reading via pandas).

If including Obspy, installing and downloading IRIS dataset more relies on internet. I can mitigate this by reducing the file size, perhaps by limiting the earthquake data to a shorter time period.

The earthquake data in CSV format is actually not too big (64KB before, now 29KB). It's more the south_america_topography_05m.txt file that's a little large (5.0MB) 🙂, but I think it's ok to just keep it in the repository for now instead of changing the workflow. Also ok with having ObsPy included in the dependency list, but that can be done in a follow-up PR (only if there's time, otherwise the PR is ok as is).

@jhtong33, I'll leave this up for a few more hours in case there's any small changes you still want to make, and will merge once you hit approve (or by the end of my work day in ~7 hours).

[weiji14]

Tutorial is now up at https://www.generic-mapping-tools.org/agu24workshop/tut04_geophysics.html 🎉