Description of the desired feature

With our Green's functions interpolators, we can predict data at any point once we have a fitted model. We could have a method BaseGridder.predict_derivative(coordinates, direction=(east, north)) to calculate (central) finite-difference derivatives on a given set of points on arbitrary directions (given by a direction vector or keywords like north and east). The method could allow for arbitrary precision (FD order) as well. This method can later be used to grid derivatives, either through a new method or an argument to BaseGridder.grid.

An example usage would be:

import verde as vd
data = vd.datasets.fetch_rio_magnetic()
coordinates = projection((data.longitude, data.latitude))
spline = vd.Spline().fit(coordinates, data.total_field_anomaly_nt)
east_deriv = spline.predict_derivative(coordinates, direction="east")
# Or using a vector
north_deriv = spline.predict_derivative(coordinates, direction=(0, 1))

A challenge would be determining the spacing for the finite differences. There are probably methods out there for estimating this (maybe Fukushima (2018)). For now, just having the spacing as an argument is more than enough. We can worry about automation later.

It would be good also to calculate second or third derivatives. This could be an argument to the method that defaults to order=1 for first-derivative. But this is not a priority and can be implemented in subsequent PRs.