GPS velocities (3-component) for the Alps

This is a compilation of 3D GPS velocities for the Alps. The horizontal velocities are reference to the Eurasian frame. All velocity components and even the position have error estimates, which is very useful and rare to find in a lot of datasets.

Original source: Sánchez et al. (2018)

Pre-processing: Source code for preparation of the original dataset for redistribution in Ensaio

import numpy as np
import pandas as pd
import pygmt

import ensaio

Download and cache the data and return the path to it on disk

fname = ensaio.fetch_alps_gps(version=1)
print(fname)

/home/runner/work/_temp/cache/ensaio/v1/alps-gps-velocity.csv.xz

Load the CSV formatted data with pandas

data = pd.read_csv(fname)
data

	station_id	longitude	latitude	height_m	velocity_east_mmyr	velocity_north_mmyr	velocity_up_mmyr	longitude_error_m	latitude_error_m	height_error_m	velocity_east_error_mmyr	velocity_north_error_mmyr	velocity_up_error_mmyr
0	ACOM	13.514900	46.547935	1774.682	0.2	1.2	1.1	0.0005	0.0009	0.001	0.1	0.1	0.1
1	AFAL	12.174517	46.527144	2284.085	-0.7	0.9	1.3	0.0009	0.0009	0.001	0.1	0.2	0.2
2	AGDE	3.466427	43.296383	65.785	-0.2	-0.2	0.1	0.0009	0.0018	0.002	0.1	0.3	0.3
3	AGNE	7.139620	45.467942	2354.600	0.0	-0.2	1.5	0.0009	0.0036	0.004	0.2	0.6	0.5
4	AIGL	3.581261	44.121398	1618.764	0.0	0.1	0.7	0.0009	0.0009	0.002	0.1	0.5	0.5
...	...	...	...	...	...	...	...	...	...	...	...	...	...
181	WLBH	7.351299	48.415171	819.069	0.0	-0.2	-2.8	0.0005	0.0009	0.001	0.1	0.2	0.2
182	WTZR	12.878911	49.144199	666.025	0.1	0.2	-0.1	0.0005	0.0005	0.001	0.1	0.1	0.1
183	ZADA	15.227590	44.113177	64.307	0.2	3.1	-0.3	0.0018	0.0036	0.004	0.2	0.4	0.4
184	ZIMM	7.465278	46.877098	956.341	-0.1	0.4	1.0	0.0005	0.0009	0.001	0.1	0.1	0.1
185	ZOUF	12.973553	46.557221	1946.508	0.1	1.0	1.3	0.0005	0.0009	0.001	0.1	0.1	0.1

186 rows × 13 columns

To plot the vectors with PyGMT, we need to convert the horizontal components into angle (azimuth) and length.

angle = np.degrees(np.arctan2(data.velocity_north_mmyr, data.velocity_east_mmyr))
length = np.hypot(data.velocity_north_mmyr, data.velocity_east_mmyr)

Now we can make a PyGMT map with the horizontal velocity vectors and vertical velocities encoded as colored points.

# West, East, South, North boundaries of the map
region = [-5, 20, 40, 55]

fig = pygmt.Figure()
with fig.subplot(
    nrows=1,
    ncols=2,
    figsize=("35c", "15c"),
    sharey="l",  # shared y-axis on the left side
    frame="WSrt",
):
    with fig.set_panel(0):
        fig.basemap(region=region, projection="M?", frame="af")
        fig.coast(area_thresh=1e4, land="#eeeeee")
        scale_factor = 2 / length.max()
        fig.plot(
            x=data.longitude,
            y=data.latitude,
            direction=[angle, length * scale_factor],
            style="v0.15c+e",
            fill="blue",
            pen="1p,blue",
        )
        # Plot a quiver caption
        fig.plot(
            x=-4,
            y=42,
            direction=[[0], [1 * scale_factor]],
            style="v0.15c+e",
            fill="blue",
            pen="1p,blue",
        )
        fig.text(
            x=-4,
            y=42.2,
            text="1 mm/yr",
            justify="BL",
            font="10p,Helvetica,blue",
        )
    with fig.set_panel(1):
        fig.basemap(region=region, projection="M?", frame="af")
        fig.coast(area_thresh=1e4, land="#eeeeee")
        pygmt.makecpt(
            cmap="polar",
            series=[data.velocity_up_mmyr.min(), data.velocity_up_mmyr.max()],
        )
        fig.plot(
            x=data.longitude,
            y=data.latitude,
            fill=data.velocity_up_mmyr,
            style="c0.2c",
            cmap=True,
            pen="0.5p,black",
        )
        fig.colorbar(
            frame='af+l"vertical velocity [mm/yr]"',
            position="jTL+w7c/0.3c+h+o1/1",
        )
fig.show()