Finding shortest distance between x,y,z point and grd surface

Q&A PyGMT Q&A
1

I want to calculate the shortest distance between an earthquake (x,y,z) and a subducting plate interface (.grd file). I.e. not just the vertical distance, but the shortest distance between a point and a curved 3D surface

Does anyone have any ideas how I could do this?

2

Maybe grdmath with one of the Xdist operator ?

3

That looks like a good starting point but as far as I can tell that will only calculate horizontal distance. So it wouldn’t consider the z value of the point (earthquake) and grd file (subducting plate)

4

I haven’t dug into it myself, but I think I’ve seen something in the gallery that might be relevant (?). The script uses sphtriangulate

5

I’m afraid we’ll need something like this https://www.mathworks.com/matlabcentral/fileexchange/33073-triangle-ray-intersection?s_tid=srchtitle

6

A bit brute force, but what if you

  1. Use grd2xyz to convert the grid surface to a bunch of XYZ points
  2. Compute the distance between each of those XYZ points to your single earthquake epicentre
  3. Sort the distances, and get the shortest distance
7

Yeah that’s what I’ve ended up doing, here’s an example if anyone else would find it useful:

import pygmt
from geopy import distance
import math

plate = pygmt.grd2xyz('plate.grd', output_type='pandas')

# this code currently loops over an obspy Catalog containing earthquake info
for event in cat:
    distances = []
    
    for index, row in plate.iterrows():
        hor_dist = distance.distance((event.origins[0].latitude, event.origins[0].longitude),
                                     (row['y'], row['x'])).km
        ver_dist = (event.origins[0].depth) - row['z']
        threeD_distance = math.sqrt(hor_dist**2 + ver_dist**2)
        distances.append(threeD_distance)

    minimum_distance = min(distances, key=abs)
8

Ah cool, exactly that! You could probably optimize things a little bit using https://geopandas.org/en/v0.11.0/docs/reference/api/geopandas.GeoDataFrame.sjoin_nearest.html and compute the distance just once, but if it’s fast enough, then don’t worry about it.