Making cross sections with project

Hi everyone,

Long time GMT user here but just trying to migrate over to PyGMT.

I want to make a cross-section of earthquakes. In GMT6 I would use something like:

awk '{print $Lon, $Lat, $Depth, $Mag}' $DATA | gmt project -C$XSectStartLon/$XSectStartLat -A$Azimuth -L$Length -Fpz -W-5/5 -Q | awk '{print $1,$2,$3}' | gmt plot -Sc -Glightgrey

It doesn’t seem like PyGMT has the project functionality, what’s the best way to approach this?

Cheers

Hi Finn,

Just seen your question as I’m also trying out PyGMT at the moment and I’ve just been tackling this exact problem!

AFAIK, with modules that haven’t been wrapped yet you have to use pygmt.clib.Session().call_module and manually fill in the options.

I’ll show you what I’ve come up with - the following function, which takes a pandas DataFrame with earthquake info (assuming first two columns are longitude and latitude) and returns the same DataFrame but with an additional column “distance”, which is the p co-ordinate output from gmt project. I’ve used the start/end (-C -E) method rather than azimuth/length, so the input variable line_coords is a list [lon_start, lat_start, lon_end, lat_end], and swath_width is the width either side of the line for input to the -W option.

def project_to_xsection(evstats, line_coords, swath_width):
    columns = evstats.columns.values.tolist()
    columns.append('distance')
    with pygmt.clib.Session() as session:
        with session.virtualfile_from_data(data = evstats) as fin:
            with pygmt.helpers.GMTTempFile() as tmp:
                session.call_module("project", f"{fin} -C{line_coords[0]}/{line_coords[1]} -E{line_coords[2]}/{line_coords[3]} -W-{swath_width}/{swath_width} -Q -Lw -Fxyzp ->{tmp.name}")
                new_coords = pd.read_csv(tmp.name, delim_whitespace = True, names = columns, dtype = 'float')
    return new_coords

Hope that gives you some idea.

Tom

1 Like

Hi Tom, thanks for that. That’s a nice work-around. I edited your code slightly to fit my dataframe format and also added an option to use azimuth and length:

def project_to_xsection_azimuth(eq_data, line_start_coords, azimuth, length, swath_width):
    tempdata = eq_data[['lon','lat','depth','mag']].copy()
    columns = tempdata.columns.values.tolist()
    columns.append('distance')
    with pygmt.clib.Session() as session:
        with session.virtualfile_from_data(data = tempdata) as fin:
            with pygmt.helpers.GMTTempFile() as tmp:
                session.call_module("project", f"{fin} -C{line_start_coords[0]}/{line_start_coords[1]} -A{azimuth} -L{length} -W-{swath_width}/{swath_width} -Q -Fxyzp ->{tmp.name}")
                projected_data = pd.read_csv(tmp.name, delim_whitespace = True, names = columns, dtype = 'float')
    return projected_data

def project_to_xsection_start_end(eq_data, line_coords, swath_width):
    tempdata = eq_data[['lon','lat','depth','mag']].copy()
    columns = tempdata.columns.values.tolist()
    columns.append('distance')

    with pygmt.clib.Session() as session:
        with session.virtualfile_from_data(data = tempdata) as fin:
            with pygmt.helpers.GMTTempFile() as tmp:
                session.call_module("project", f"{fin} -C{line_coords[0]}/{line_coords[1]} -E{line_coords[2]}/{line_coords[3]} -W-{swath_width}/{swath_width} -Q -Fxyzp ->{tmp.name}")
                projected_data = pd.read_csv(tmp.name, delim_whitespace = True, names = columns, dtype = 'float')
    return projected_data

Cheers,
Finn

1 Like

Nice!