GeoMakie.jl: Stereographic projection at pole did not work

Hi all,

I started to code in Julia half a year ago and I am quite excited about using GeoMakie. I work on ice sheet dynamics so I mostly want to visualize the poles. Therefore I want to use a stereographic projection and as a first example I tried to adapt the keyword argument in a GeoMakie example of the documentation (Examples · GeoMakie.jl). The problem I face is that I obtain a blank plot:

Here goes the code:

using GeoMakie, CairoMakie

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

fig = Figure()
ga = GeoAxis(
    fig[1, 1],
    dest="+proj=sterea +lat_0=90",    # Only lined that was changed compared to example!
    lonlims = automatic,
    coastlines = true,
    title = "Stereographic projection"
)

sp = surface!(ga, lons, lats, field; shading = false, colormap = :rainbow_bgyrm_35_85_c69_n256)
cb = Colorbar(fig[1, 2], sp)

fig

I’d be very thankful if anyone can help!

Just wondering whether that’s the correct spelling… ?

Should exist:

https://proj.org/operations/projections/sterea.html

Maybe @asinghvi17 can help

:slight_smile: well you learn something every day

If you want to consider an alternative

using GMT
lons = -180:180;
lats = 0:90;      # No need for the whole Earth
field = [exp(cosd(l)) + 3(y/90) for y in lats, l in lons];

G = mat2grid(field, x=lons, y=lats, proj4="+proj=longlat");
imshow(G, proj="+proj=sterea +lat_0=90", coast=true)

Thanks a lot for the neat answer! On the longer term I’d really like to rely as much as possible on native julia tools, so I’d still be very interested if anyone has a solution in GeoMakie :slight_smile:

That is an illusion. No plotting package is Julia native.

Well it does make a difference if the primitives are julia objects because that’s what people mainly hack around with. Of course the rendering itself is not Julia.

1 Like

I suspect that with that definition then all (most?) plotting packages are Julia natives.

They differ in the depth that their Julia stack reaches

Then I guess my wish would be better formulated by saying: I want to use (Geo)Makie because it presents a large julia stack and it would be more coherent with the rest of my plotting routines!

Maybe because of those gotchas=> Redirecting to https://geo.makie.org/stable or => Redirecting to https://geo.makie.org/stable ?

I see no reply so probably there is some other cause that makes the plot not showing. Just wanted to add that I like both of those plotting packages (GMT and GeoMakie).

1 Like

Following @j_u’s suggestion on one of the gotchas. Playing with the xlimits! and ylimits! make the plot appear.


using GeoMakie, CairoMakie

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

fig = Figure()
ga = GeoAxis(
    fig[1, 1],
    dest="+proj=sterea +lat_0=90",    # Only lined that was changed compared to example!
    lonlims = automatic,
    coastlines = true,
    title = "Stereographic projection"
)

xlims!(ga, -90, 90)
ylims!(ga, -50, 50)
sp = surface!(ga, lons, lats, field; shading = false, colormap = :rainbow_bgyrm_35_85_c69_n256)
cb = Colorbar(fig[1, 2], sp)
fig

Produces:

1 Like

Muchas gracias Argel, it solved it!

Here goes my code (this time for Antarctica), where the aspect ratio is fixed and the limits are adapted. Hope it can be useful for others.

using GeoMakie, CairoMakie

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

fig = Figure()
ga = GeoAxis(
    fig[1, 1],
    aspect = AxisAspect(1),
    dest="+proj=sterea +lat_0=-90",
    lonlims = automatic,
    coastlines = true,
    title = "Stereographic projection"
)

xlims!(ga, -180, 180)
ylims!(ga, -90, -60)
sp = surface!(ga, lons, lats, field; shading = false, colormap = :rainbow_bgyrm_35_85_c69_n256)
cb = Colorbar(fig[1, 2], sp)
fig

Should display the following picture:

2 Likes

Great to see it working on your sides @aramirezreyes and @JanJereczek and sorry for providing only a very general advice. Really happy that GeoMakie as I understand received a significant refresh recently. Hope the development efforts will continue.

Hey, was on vacation for a bit so just saw this. Thanks @aramirezreyes for the solution!

Yeah, in general this is an issue in GeoMakie - if the limits of the plot are greater than the intrinsic limits of the projection, the whole thing goes blank, since the transformed limits become infinity. It’s hard to define a general solution because not all projections have well-defined rectangular boundingboxes.

You will see that the surface is not clipped by the edge of the Axis. My current thought is to just add a white polygon on top to simulate clipping until we get a proper nonlinear clip implemented, but not totally sure how that should work still.