Plotting vectors on oblique mercator map

[ledeczi]

Hi all,
I am plotting a map using an oblique mercator projection. I’d like to add vectors with earthquake latitude, longitude, and nodal plane strike. However, when I plot these as vectors on top of the map, they don’t plot in the oblique mercator projection of the basemap—they plot using the cartesian coordinates of the figure, so the azimuths are wrong. Here is a photo of the comparison between the same azimuths plotted on the two projections. The earthquakes should strike roughly E-W the entire western segment, so they should not also be striking that same direction in the oblique mercator map. This is even after I specify the projection when plotting these vectors.

See code:

fig = pygmt.Figure()
fig.coast(projection="OC147/-3/-10.7/145.7/6c", region="141/-4.0/155/-3.0+r",
    frame="afg",
    land="gray",
    shorelines="1/thin",
   water="lightblue",
)

fig.plot(region="141/-4.0/155/-3.0+r",
    projection="OC147/-3/-10.7/145.7/6c",
    data=vectors,
    style="v0c",
    pen="0.1p",
)
fig.show(width=1000, dpi=400)

#vectors are in format [latitude, longitude, azimuth, length]

[yvonnefroehlich]

Hello @ledeczi,

Welcome to the GMT forum !

If you want to plot geographic vector you need to use = together with the style parameter. We have a gallery example at Cartesian, circular, and geographic vectors — PyGMT and a tutorial at Plotting vectors — PyGMT, which show how to plot geographic vectors.

If it is possible for you please also provide you input data “vectors”, so people can reproduce your issue.

[ledeczi]

Hi Yvonne, Thanks so much for your reply!

Here is an example of 5 of my vectors:

       [144.28 ,  -3.407, 359.8  ,   0.4  ],
       [144.43 ,  -3.412, 361.4  ,   0.4  ],
       [144.713,  -3.418, 363.   ,   0.4  ],
       [144.757,  -3.431, 362.1  ,   0.4  ],
       [144.949,  -3.408, 361.3  ,   0.4  ],

I saw the geographic coordinates, but saw that it requires a start and end point. I can do this calculation but it would be easier and more accurate to input my azimuths rather than calculating an endpoint. Is that the easiest solution or is there a way to use my azimuths as-is? If not, that’s ok!

[yvonnefroehlich]

Yes, it’s also possible to use the azimuth and length of geographic vectors for plotting. But please note, that the length is the geographic length (by default in kilometers).

import pygmt as gmt
import numpy as np


lons = [-80, -60, 0, -40, -20]
lats = [0] * len(lons)
vector_dirs = [0, 10, 20, 30, 40]
vector_lens = np.array([1, 2, 3, 4, 5]) *1000 

# lon_degE, lat_degN, dir_cw_from_N, len_km
vectors = [
    [-80, 0,  0, 1000], 
    [-60, 0, 10, 2000],
    [  0, 0, 20, 3000],
    [-40, 0, 30, 4000],
    [-20, 0, 40, 5000],
]

fig = gmt.Figure()

fig.plot(
    region="d",
    projection="N10c",
    frame="afg",
    data=vectors,
    style="=0.2c+e+h0.1c",
    pen="0.3p,red3",
    fill="red3",
)

fig.show()