Subscribe to our weekly newsletter

How to plot great circle path through your region using PyGMT

Utpal Kumar   3 minute read      

In this article, we will learn how to visualize the great circle paths that traverse a designated region of interest using PyGMT. Understanding these paths is crucial for rigorous seismic tomography investigations, as they can provide insights into the subsurface structure below the select region of interest.

Introduction

Seismic tomography is an indispensable tool in the study of the Earth’s internal structure. It enables researchers to develop high-resolution models of subsurface geology, thereby facilitating our understanding of tectonic processes, resource exploration, and even earthquake hazard assessment. However, to fully grasp the complexity of the Earth’s interior, one must not only focus on the source and receiver locations but also the paths that seismic waves travel between them, often termed as “great circle paths.”

This article aims to discuss a Python-based approach to visualize great circle paths that seismic waves follow, focusing on a designated region of interest. We will utilize PyGMT (Python Generic Mapping Tools), a powerful library for creating high-quality geographical maps and plots, along with other Python libraries such as NumPy, Pandas, and Shapely for data manipulation and geometrical operations.

By the end of this article, you will have learned how to use Python code to:

  • Define a specific geographic region of interest.
  • Generate and plot the great circle paths between earthquake events and seismic stations.
  • Filter the relevant stations based on their location and the paths that traverse the specified region.

Understanding these paths is not just an academic exercise; it is crucial for rigorous seismic tomography investigations. The generated paths can provide valuable insights into the subsurface structure below the selected region, thereby making your seismic analyses more robust and accurate.

Install PyGmt

Using python env

python -m venv geoviz
source geoviz/bin/activate
pip install pygmt

For more details, visit the article A Quick Overview on Geospatial Data Visualization using PyGMT.

Plot great-circle path traversing the region of interest

This Python script serves as a specialized tool for visualizing seismic events and their associated seismic stations within a defined geographical region. Utilizing the PyGMT library for geospatial plotting, along with data manipulation libraries like NumPy and Pandas, the script reads earthquake event and seismic station data, filters them based on their location, and plots them on a high-resolution map. It specifically highlights the great circle paths between the earthquake source and the seismic stations, aiding in seismic tomography studies by offering insights into the subsurface structure beneath the region of interest. Through the use of geometrical libraries like Shapely and geodetic calculations from Pyproj, the script efficiently determines whether these paths traverse the specified region, thereby streamlining the process of selecting relevant data for further analysis.

import pygmt
import numpy as np
import pandas as pd
import yaml, glob, os, sys
from shapely.geometry.polygon import Polygon
import pyproj
from shapely.geometry import Point



def get_region_polygon(
    # Box size
    lon_left   = -128,      # possible range: -180, 180 deg
    lon_right  = -96,      # possible range: -180, 180 deg
    lat_bottom = 27,      # possible range: -90, 90 deg
    lat_top    =  52,      # possible range: -90, 90 deg
    offset = 20, # offset in degrees from the box limits
):
    lon_left  = lon_left - offset
    lon_right = lon_right + offset
    lat_bottom = lat_bottom - offset
    lat_top = lat_top + offset

    box_lims = [[lon_left,lat_bottom], [lon_right,lat_bottom], [lon_right,lat_top], [lon_left,lat_top], [lon_left,lat_bottom]]
    box_maxdim = max(np.abs(lon_right-lon_left),np.abs(lat_top-lat_bottom))
    lims_array = np.array(box_lims)
    boxclon, boxclat = np.mean(lims_array[:, 0]),np.mean(lims_array[:, 1])
    box_polygon = Polygon(box_lims)

    return box_polygon

def is_in_domain(lon_points, lat_points, box_polygon):
    
    for lat, lon in zip(lat_points, lon_points):
        # print(lon, lat)
        if box_polygon.contains(Point(lon,lat)):
            return True
    return False


def main():     

    ## Inversion domain
    lon_left = -128      # possible range: -180, 180 deg
    lon_right = -96      # possible range: -180, 180 deg
    lat_bottom = 27      # possible range: -90, 90 deg
    lat_top =  52      # possible range: -90, 90 deg
    
    ## Event info
    evbase = 'C201210240045A' # (D12O0TUA)
    evlon = -85.30
    evlat = 10.09
    evdep = 17.0
    evmag = 6.0


    ## Define polygon
    box_polygon = get_region_polygon(offset = 5)
    print("box_polygon.bounds: ",box_polygon.bounds)
    geod=pyproj.Geod(ellps="WGS84")

    out_image_station = f"{evbase}.png"

    ## plot stations
    fig = pygmt.Figure()
    projection = "W-110.5885/12c"
    fig.basemap(region='g', projection=projection, frame=["afg", f"+t{evbase}"])
    fig.coast(
        land="lightgrey",
        water="white",
        shorelines="0.1p",
        frame="WSNE",
        resolution='h',
        area_thresh=10000
    )
    
    if is_in_domain([evlon], [evlat], box_polygon):
        print(f"--> Skipping {evbase} because it is in domain")
    
    dff_event = pd.read_csv('event_station_info_D12O0TUA.txt', sep='\s+', header=None, names=['evname', 'stn', 'slon', 'slat'])
    # print(dff_event)

    assert len(dff_event) > 0, "No stations found in event_station_info_D12O0TUA.txt"

    fig.plot(x=evlon, y=evlat, style="c0.3c", color="red", pen="black")

    for stlat, stlon, sname in zip(dff_event.slat, dff_event.slon, dff_event.stn):
        if is_in_domain([stlon], [stlat], box_polygon):
            continue
        line_arc=geod.inv_intermediate(evlon,evlat,stlon,stlat,npts=300)
        lon_points=np.array(line_arc.lons)
        lat_points=np.array(line_arc.lats)
        if not is_in_domain(lon_points, lat_points, box_polygon):
            continue
        fig.plot(x=lon_points, y=lat_points, pen="0.5p,black")
        fig.plot(x=stlon, y=stlat, style="i0.1c", color="blue", pen="black")
    

        rectangle = [box_polygon.bounds]
        fig.plot(data=rectangle, style="r+s", pen="2p,red")

    print('----> Saving map... {}'.format(out_image_station))
    fig.savefig(out_image_station, crop=True, dpi=600)

if __name__ == "__main__":
    main()

Download the event_station_info_D12O0TUA.txt from here

Great circle paths for event C201210240045A crossing the region of interest
Great circle paths for event C201210240045A crossing the region of interest

Disclaimer of liability

The information provided by the Earth Inversion is made available for educational purposes only.

Whilst we endeavor to keep the information up-to-date and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.

UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.

Leave a comment