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

Q&A PyGMT Q&A
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)