very slow geoid calculation

[alex-s-gardner]

I recently issued a pull request to add geoid function to Geodesy.jl because the geoid conversion in Proj.jl is very slow for large numbers of points [~>500] relative to reading the geoid from a local file. I'm not sure if Geodesy.jl is the right home for such a function or if we just need to improve Proj.jl.

Here is a figure and the corresponding code that demonstrates just how slow using Proj is for geoid to ellipsoid conversions:

# retrieve EGM2008 geoid heights ["EPSG:3855"] above WGS84 ellipsoid ["EPSG:4979"] 
# for `n random points

using GeoArrays, Proj

# file can be downloaded here: 
# https://s3-eu-west-1.amazonaws.com/download.agisoft.com/gtg/us_nga_egm2008_1.tif
geoid_file = "/Users/gardnera/data/geoids/us_nga_egm2008_1.tif"
n = [2, 10, 100, 5000, 10000];
time_local = zeros(size(n));
time_proj = zeros(size(n));

height1 = [];
height2 = [];
for i = eachindex(n)
    lat = (rand(n[i]) .- 0.5) .* 180;
    lon = (rand(n[i]) .- 0.5) .* 360;

    ## Approach 1: read directly from 1/60 degree geotiff
    begin
        t = time();
        # read [lazzily] in EGM2008 geoid as a GeoArray with 1 minute resolution
        
        egm2008 = GeoArrays.read(geoid_file);

        # find indices of points into cropped dem 
        xy = hcat(lon,lat);
        xy = [copy(r) for r in eachrow(xy)];
        ind = indices.(Ref(egm2008), xy);
        ind = [CartesianIndex((ij[1], ij[2])) for ij in ind];

        # extract values
        height1 = egm2008[ind];
        time_local[i] = time()-t;
    end;

    ## Approach 2: use Proj.jl
    begin
        t = time();
        # enamble network for geoid conversions
        Proj.enable_network!();

        # build transformation 
        trans = Proj.Transformation("EPSG:4326", "EPSG:3855", always_xy=true);

        # project points
        data = trans.(lon,lat,0);

        height2 = -getindex.(data, 3);
        time_proj[i] = time()-t;
    end;
end


using CairoMakie
# Check that both return same resolution
f = Figure()
ax = Axis(f[1, 1],
    title = "validation",
    xlabel = "local",
    ylabel = "Proj.jl"
)
scatter!(ax, height1, height2)

# plot time between appraoches
ax = Axis(f[1, 2],
    title = "timing",
    xlabel = "number of points",
    ylabel = "time [s]"
)
lines!(ax, n, time_local, label = "local")
lines!(ax, n, time_proj, label = "Proj.jl")
axislegend(ax)
f