API for geoplots with GeoMakie.jl and docs

[Datseris]

This is the initial sketch of how I, a scientist working daily with climate data, would like to use a geospatial plotting package.

The interface is showcased at examples/new_api.jl. All source code is in the src/api.jl file. Now GeoMakie.jl uses and reexports Proj4. I never understood why this wasn't the case from the get-go, as GeoMakie.jl only makes sense to be used with Proj4.jl.

Todos:

• Incorporate GeoAxis from Template for GeoAxis  #92
• Solve the GLMakie issue: even hovering over the screen yields errors, while any interactivity is broken
• There is a problem with contourf: it is incorrectly zoomed.
• There is a problem with surface: it's interpolation is incorrect and leads to artifacts.
• Finish writing docs
• Make docs runnable and presentable locally using Documenter
• Correct z-ordering of coastlines and axis grid lines so that they are guaranteed above other plot elements
• add an actual docstring for geo2basic (I don't know what it does)
• Update README

What will be left for future PRs:

• Automatic axis ticks (for better support for local plots)
• Function for "antimeridian cutting"
• Coastlines plot uses antimeridian cutting
• Add inverse transform from Pro4 (@visr )
• Add actual, proper tests

[Datseris]

output of geosurface:

output of geoscatter:

[Datseris]

@SimonDanisch can you please help me? I've created the GeoAxis, but it seems like the projection isn't actually used. I now get images like

which ar enot in this winkel tripel projection... :(

[Datseris]

I added gridlines but the projection is still not working, I don't know why. I'm using exactly the same code what the fuck. @lazarusA maybe you see the mistake?

[Datseris]

ok, it is fixed. Please teach me on how to set the axis limits because something tells me the code from lazaro is an inefficient way to do this.

Here is the status quo:

I request a review before starting working on the docs @SimonDanisch

[Datseris]

GLMakie.jl interactivity does not work here. I canot zoom in and I get errors even when hovering over the figure.

[lazarusA]

As for the axis limits, yes we only need 4 points, actually I think is just 2 points, the lower-left corner and the upper-right corner (the output from FRect2D). My issue here from the start was that these limits are different for each projection. Namely, what is the appropriate combination of (Lon,Lat)[left-corner] and (Lon-Lat)[right-corner] that will give you the minima and maxima when projected in the 2D axis. These projection information should be somewhere (it should be known), but I haven't found a reference. Hence, my approach is very naive and inefficient.

[Datseris]

Yeap you are right and I was wrong. You need much more points for a generic projection, as you can't know the max/min lon/lat just from the "corners" !

[codecov-commenter]

Codecov Report

Merging #93 (1593567) into master (276b90d) will decrease coverage by 15.13%.
The diff coverage is 0.00%.

@@             Coverage Diff             @@
##           master      #93       +/-   ##
===========================================
- Coverage   35.13%   20.00%   -15.14%     
===========================================
  Files           2        3        +1     
  Lines          37       65       +28     
===========================================
  Hits           13       13               
- Misses         24       52       +28     

Impacted Files	Coverage Δ
src/api.jl	0.00% <0.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 276b90d...1593567. Read the comment docs.

[Datseris]

People, there is a problem. This doesn't work for contourf!:

lons = -180:180
lats = -90:90
field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

# Surface example
fig = Figure()
ax = GeoAxis(fig[1,1])
el = surface!(ax, lons, lats, field)
display(fig)

The internal logic of surface! and contourf! seem to be different w.r.t. how they handle axis transforms?

[Datseris]

Actually, surface! is also wrong. It incorrectly plots data. Using a physical dataset that I am very familiar with it, I've tested it, and it shows very wrong stuff. Using a scatter plot, this is how the physical data look like:

Now, plotting the same data as a surface plot, I get:

there are many horizontal artifacts that are incorrect. :( I think we need to fix contourf but I don't know how...

[SimonDanisch]

Hm contourf seems to have a zooming issue:

I'm not really sure why...I'd suspect the boundingbox, but that's the same as for surface

[lazarusA]

Suggested change

# Plot coastlines

[lazarusA]

Suggested change

if coastlines

[lazarusA]

Suggested change

coastplot = lines!(ax, GeoMakie.coastlines(); color = :black, overdraw = true, coastkwargs...)

[lazarusA]

Suggested change

translate!(coastplot, 0, 0, 99) # ensure they are on top of other plotted elements

[lazarusA]

Suggested change

end

[lazarusA]

@Datseris Since you agree with the coastlines arg, I also added suggestions. You can batch the suggestions you agree on and then commit them here directly. No need for additional commit/push from local.

[SimonDanisch]

@Datseris do you want to update this PR with the suggestions? Would be nice to merge it this week finally :)

[Datseris]

@SimonDanisch do not forget that there are some todos for you here lingering for a while now. Specifically the incorrect "coastlines always at the top of the plot", the "infinite error messages when hovering over window at GLMakie" and the "countrf plot is plotted zoomed in".

[Datseris]

and someone has to solve:

add an actual docstring for geo2basic (I don't know what it does)

[gaelforget]

Of course it's better to first get it working, but I don't want to forget about the local/regional usecases.

Agreed. Might add the need for stereographic polar projections as a common use case -- anyone working on the Arctic or Antarctic regions will want those, for good reasons.

Here is an example of what I used to do in a Matlab using https://www.eoas.ubc.ca/~rich/map.html for global ocean simulations. Two stereo projection + one map where I cirshifted the dateline such that 20E is the right most limit. Ideally I would like to be able to do this with Makie.jl & coastlines.

[Datseris]

This plot is already possible to my knowledge, provided you know the projection's name for Proj4. What is NOT possible at the moment is the longitude and latitude labels. At the moment their implementation is generic and I have a suspicion they wont work well for specialized projections.

[Datseris]

Oh, what is also not possible is to make the continents gray. At least I don't know how. I only have seen how to plot the coastlines but not fill them in. Ah, also, coastline plotting does not work in most projections, only those that are somewhat rectangular and go from -180 to 180 degrees.

[gaelforget]

What is NOT possible at the moment is the longitude and latitude labels.

Indeed. Highlighting this was one of the reasons to post this image. The circular shape of the domain edge, which I dont know how to deal with either, is another potential hurdle.

At the moment their implementation is generic and I have a suspicion they wont work well for specialized projections.

any projection is specialized I feel -- including the use of lon & lat as x & y

[gaelforget]

Oh, what is also not possible is to make the continents gray. At least I don't know how. I only have seen how to plot the coastlines but not fill them in.

Not sure how m_map does it (or the lakes showing up in black)

Ah, also, coastline plotting does not work in most projections, only those that are somewhat rectangular and go from -180 to 180 degrees.

As @visr pointed out the only reason coastlines work with the -180 to 180 maps is the polygons cuts done per the GeoJSON specs. Don't know how to generalize this either for an arbitrary domain shape (e.g. regional, polar, or rotated regional)

[lazarusA]

Ok. This works (almost all the time :D)

using GeoMakie, CairoMakie
using Proj4, GeometryBasics, GeoDatasets

lonl, latl, data = GeoDatasets.landseamask(; resolution = 'c', grid = 10)
cmap = cgrad(:Greys_4, 3, categorical = true, rev = false)

fig = Figure(resolution = (1200, 800), fontsize = 22)
ax = Axis(fig[1, 1], aspect = DataAspect(),
    title = "Projection: airy")
# all input data coordinates are projected using this function
trans = Proj4.Transformation("+proj=longlat +datum=WGS84", "+proj=airy +lat_0=0 +lon_0=0",
    always_xy = true)
# "+proj=laea +lat_0=-45 +lon_0=-90"
# "+proj=stere +lat_0=10 +lon_0=0"
# see https://proj.org/operations/projections/index.html for more options 
ptrans = Makie.PointTrans{2}(trans)
ax.scene.transformation.transform_func[] = ptrans
# add some limits, still it needs to be manual
lats = -90:30.0:90 |> collect
lons = -180:30.0:180 |> collect
# avoid PROJ wrapping 180 to -180
lons[end] = prevfloat(lons[end])
points = [Point2f0(lon, lat) for lon in lons, lat in lats]
xytrans = Makie.apply_transform(ptrans, points)
function checkInfs(xy)
    xy != [Inf, Inf] && xy != [-Inf, -Inf] && xy[1] != Inf && xy[2] != Inf && xy[1] != -Inf && xy[2] != -Inf
end
xytransClean = [xytrans[i] for i in 1:prod(size(xytrans)) if checkInfs(xytrans[i])]
rectLimits = FRect2D(xytransClean)
# now the plot
pltobj = heatmap!(ax, lonl, latl, data, colormap = cmap)
lines!(ax, GeoMakie.coastlines(), color = :black)
limits!(ax, rectLimits)
# optional
lonlines = [Point2f0(j, i) for i = -90:90, j in lons]
latlines = [Point2f0(j, i) for j = -180:180, i in lats]
[lines!(ax, lonlines[:, i], color = (:black, 0.25),
    linestyle = :dot, overdraw = true) for i = 1:size(lonlines)[2]]
[lines!(ax, latlines[:, i], color = (:black, 0.25), linestyle = :dot,
    overdraw = true) for i = 1:size(latlines)[2]]

hidedecorations!(ax, ticks = false, label = false)
hidespines!(ax)
save("airy.png", fig)
fig

And for the stereo projection heatmap gives problems, but surface works[kinda] (it will be a matter of setting proper limits for each region afterwards):

fig = Figure(resolution = (1200, 800), fontsize = 22)
ax = Axis(fig[1, 1], aspect = DataAspect(),
    title = "Projection: stere")
# all input data coordinates are projected using this function
trans = Proj4.Transformation("+proj=longlat +datum=WGS84", "+proj=stere",
    always_xy = true)
# "+proj=laea +lat_0=-45 +lon_0=-90"
# "+proj=stere +lat_0=10 +lon_0=0"
# see https://proj.org/operations/projections/index.html for more options 
ptrans = Makie.PointTrans{2}(trans)
ax.scene.transformation.transform_func[] = ptrans
# add some limits, still it needs to be manual
lats = -90:30.0:90 |> collect
lons = -180:30.0:180 |> collect
# avoid PROJ wrapping 180 to -180
lons[end] = prevfloat(lons[end])
points = [Point2f0(lon, lat) for lon in lons, lat in lats]
xytrans = Makie.apply_transform(ptrans, points)
function checkInfs(xy)
    xy != [Inf, Inf] && xy != [-Inf, -Inf] && xy[1] != Inf && xy[2] != Inf && xy[1] != -Inf && xy[2] != -Inf
end
xytransClean = [xytrans[i] for i in 1:prod(size(xytrans)) if checkInfs(xytrans[i])]
rectLimits = FRect2D(xytransClean)
# now the plot
lats = -90:1.0:90 |> collect
lons = -180:1:180 |> collect
field = [exp(cosd(l)) + 3(y / 90) for l in lons, y in lats]
hm1 = surface!(ax, lons, lats, field, shading = false, overdraw = true)
#pltobj = heatmap!(ax, lonl, latl, data, colormap = cmap)
lines!(ax, GeoMakie.coastlines(), color = :black)
limits!(ax, rectLimits)
# optional
lats = -90:30.0:90 |> collect
lons = -180:30.0:180 |> collect

lonlines = [Point2f0(j, i) for i = -90:90, j in lons]
latlines = [Point2f0(j, i) for j = -180:180, i in lats]
[lines!(ax, lonlines[:, i], color = (:black, 0.25),
    linestyle = :dot, overdraw = true) for i = 1:size(lonlines)[2]]
[lines!(ax, latlines[:, i], color = (:black, 0.25), linestyle = :dot,
    overdraw = true) for i = 1:size(latlines)[2]]

hidedecorations!(ax, ticks = false, label = false)
hidespines!(ax)
save("stere.png", fig)
fig

[Datseris]

the coastlines and the "land sea mask" are different though, there are many places where they do not match well.

[Datseris]

In the second code block, the only thing that is different from current source code is the chcking for infinite, right?

[lazarusA]

Yes. They are two different data sets. So is to be expected. And yes, I added an extra step to check for Infs.

[SimonDanisch]

@SimonDanisch do not forget that there are some todos for you here lingering for a while now. Specifically the incorrect "coastlines always at the top of the plot", the "infinite error messages when hovering over window at GLMakie" and the "countrf plot is plotted zoomed in".

So, I just implemented some very basic depth sorting for GLMakie, but I can't find a runnable example reproducing the bug you have... @Datseris could you post one?

infinite error messages when hovering over window at GLMakie"

That's fixed on master

[SimonDanisch]

Contourf resets the limits, because it uses tight_limits:

fig = Figure()
ax = Axis(fig[1,1], limits = ((0, 2), (0, 2)))
contourf!(ax, -1..1, -1..1, rand(4, 4)) # resets limits to -1..1, ignores set limits from axis
fig

I'm talking with @jkrumbiegel how to resolve this...
Ideally, we teach the Axis to correctly calculate the limits with the point transform, and fix the tight_limit bug

[SimonDanisch]

I think I fixed the limits, so that one doesn't need to calculate them in GeoMakie... Although I remember faintly, that the fix I used may make some problems with other examples...

Well, I'll see if the tests pass...

[Datseris]

ok, so I just need to test if they layering of lines is correct. Okay will do that asap

[Datseris]

continued in #95