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

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?

Maybe grdmath with one of the Xdist operator ?

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)

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

I’m afraid we’ll need something like this

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

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(([0].latitude,[0].longitude),
                                     (row['y'], row['x'])).km
        ver_dist = ([0].depth) - row['z']
        threeD_distance = math.sqrt(hor_dist**2 + ver_dist**2)

    minimum_distance = min(distances, key=abs)

Ah cool, exactly that! You could probably optimize things a little bit using and compute the distance just once, but if it’s fast enough, then don’t worry about it.