Use PyGMT for GLOBK results

Q&A PyGMT Q&A
1

Hello,
I’m a french student doing an internship on GNSS data, and I have to use Gamit/Globk to process it. I’m able to generate a map with displacement vectors, but currently, this map(made by sh_plotvel function) is all white with only the coastlines in black. I would like to add a topographic background.

I can create a map with a topographic background separately, but I can’t figure out how to incorporate this background into “sh_plotvel” in the code below. If anyone has any ideas on how to do this, I would be very grateful.

Thank you.

# Create the map using PyGMT
fig = pygmt.Figure()
fig.coast(
    region=[0, 9, 42, 47],  # Region for France
    projection="H15c",   # Lambert 93 projection
    borders=1,
    land="lightgray",
    water="lightblue",
    area_thresh=10000,
    shorelines="0.5p,black",
    frame=["a", "+tStation positions"],
)

# Add the relief
grid = pygmt.datasets.load_earth_relief(
    resolution="30s", region=[0, 9, 42, 47], registration="gridline"
)
fig.grdimage(grid=grid, frame="a", projection="H15c", cmap="oleron", transparency=20)

# Add sh_plotvel
map_name = input("Please enter the map name: ")
sh_plotvel_args = [
    'sh_plotvel',
    '-ps', map_name,
    '-f', 'globk_vel.org',
    '-factor', '1',
    '-arrow_value', '20',
    '-R0/9/42/47',
]
subprocess.run(sh_plotvel_args, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)

# Add the legend, scale, and compass rose
with pygmt.config():
    fig.basemap(rose="jT+w2.5c+lO,E,S,N+o1c/0.3c", map_scale="jBL+w300k+o0.5c/0.5c+f")

for position in ("T"):
    fig.text(
        text="Lambert 93 Projection",
        x=6.3,
        y=41.9,
        font="10p,Helvetica-Bold,black",
    )

for position in ("T"):
    fig.text(
        text="Km",
        x=-2.8,
        y=42.3,
        font="10p,Helvetica-Bold,black",
    )

# Show the legend
fig.legend(position="JMR+o1c/0c", box=True)

# Show the map
fig.show()
print('\nsh_org2vel output')
sh_org2vel_args = ['sh_org2vel', '-file', 'globk_comb.org']
if t == 1:
    subprocess.run(sh_org2vel_args, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)

print('\nStarting GLobk')
globk_args = ['globk', '6', 'globk_comb.log', 'globk_comb.prt', ''+nom_experience+'.gdl', 'globk.cmd']
if t == 1:
    subprocess.run(globk_args, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)

# Execute command to convert PS files to JPEG images
print('\nConverting PS files to JPEG images')
os.chdir(vsoln)
psconvert_args = ['gmt psconvert *.ps -Tj']
subprocess.run(psconvert_args, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, shell=True)

print('\nProgram finished')
2

Hi,

this thing sh_plotvel seems using GMT internally for plotting according to its documentation. I could not check the code as the software supplier requires registration for downloading.

Pygmt uses a mechanism called gmt session. This is to ensure pygmt’s plotting calls’ output lands into the same figure. sh_plotvel is totally unaware of this. The result is pygmt calls do its own figure and sh_plotvel does its own. Basically you should get two figures out of your script, one from pygmt calls and another from sh_plotvel.

To combine sh_plotvel with pygmt or command-line gmt plotting, you likely need to check its code, extract and modify its plotting calls it to use with command-line gmt plotting, or to port the plotting calls to pygmt to use it with pygmt plotting.

While possible, it may not be super easy