Add EQL gridder in geodetic coordinates

[santisoler]

Description of the desired feature

Would be nice to have an equivalent layer gridder that works in geodetic coordinates.
This means that the data points and the location of the point sources should be defined in geodetic coordinates.
The Green's function should still be equal to 1 over the Euclidean distance between each data point and each source.
After #154 is closed, we will be able to easily compute this distance between two given points defined in geodetic coordinates.
The new EQL gridder should also take the reference ellipsoid, which should be a instance of boule.Ellipsoid. By default, we should set it to the WGS84 ellipsoid, for which we need to add Boule as a dependency of Harmonica.
Would be better to wait until #165 is fixed before start working on this issue: the idea is to have a common base class for every EQL gridder.

After having a geodetic EQL, we could remove the gallery example for the spherical EQL.
Using the geodetic EQL will be more straight forward for users (they don't need to perform coordinate conversions), but we should keep the spherical EQL on the API just in case we want it in the future.

Are you willing to help implement and maintain this feature? Yes