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
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