Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

API for geoplots with GeoMakie.jl and docs #93

Closed
wants to merge 35 commits into from

Conversation

Datseris
Copy link
Collaborator

@Datseris Datseris commented Oct 29, 2021

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
Copy link
Collaborator Author

output of geosurface:
image

output of geoscatter:
image

@Datseris Datseris mentioned this pull request Oct 29, 2021
@Datseris
Copy link
Collaborator Author

@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
image

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

@Datseris
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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:
image

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

@Datseris
Copy link
Collaborator Author

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

@lazarusA
Copy link
Contributor

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
Copy link
Collaborator Author

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
Copy link

codecov-commenter commented Oct 30, 2021

Codecov Report

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

Impacted file tree graph

@@             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
Copy link
Collaborator Author

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)

image

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

@Datseris
Copy link
Collaborator Author

Datseris commented Oct 31, 2021

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:
image

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

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

@SimonDanisch
Copy link
Member

Hm contourf seems to have a zooming issue:
image
I'm not really sure why...I'd suspect the boundingbox, but that's the same as for surface

src/geoaxis.jl Outdated Show resolved Hide resolved
src/geoaxis.jl Outdated Show resolved Hide resolved
src/geoaxis.jl Outdated Show resolved Hide resolved
rectLimits = FRect2D(Makie.apply_transform(ptrans, points))
limits!(ax, rectLimits)

# Plot coastlines
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# Plot coastlines

limits!(ax, rectLimits)

# Plot coastlines
if coastlines
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if coastlines


# Plot coastlines
if coastlines
coastplot = lines!(ax, GeoMakie.coastlines(); color = :black, overdraw = true, coastkwargs...)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

# Plot coastlines
if coastlines
coastplot = lines!(ax, GeoMakie.coastlines(); color = :black, overdraw = true, coastkwargs...)
translate!(coastplot, 0, 0, 99) # ensure they are on top of other plotted elements
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

if coastlines
coastplot = lines!(ax, GeoMakie.coastlines(); color = :black, overdraw = true, coastkwargs...)
translate!(coastplot, 0, 0, 99) # ensure they are on top of other plotted elements
end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
end

Copy link
Contributor

@lazarusA lazarusA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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
Copy link
Member

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

Datseris and others added 2 commits December 15, 2021 22:20
Co-authored-by: Lazaro Alonso <lazarus.alon@gmail.com>
@Datseris
Copy link
Collaborator Author

@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
Copy link
Collaborator Author

and someone has to solve:

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

@gaelforget
Copy link

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.

Screen Shot 2021-12-15 at 4 43 59 PM

@Datseris
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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
Copy link

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
Copy link

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
Copy link
Contributor

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

airy

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

stere

@Datseris
Copy link
Collaborator Author

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

@Datseris
Copy link
Collaborator Author

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

@lazarusA
Copy link
Contributor

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

@SimonDanisch
Copy link
Member

@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
Copy link
Member

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
Copy link
Member

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...
image
Well, I'll see if the tests pass...

@Datseris
Copy link
Collaborator Author

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

@Datseris
Copy link
Collaborator Author

continued in #95

@Datseris Datseris closed this Dec 22, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants