Using GeoMakie

[Datseris]

Here is a code example I use to animate my fields over the earth (and also do a single plot of their time averages). In the following assume that A is a 3D matrix with dimensions lon x lat x t , with lon lat in degrees and t in months. Because my lon data are structured starting in 0 and ending in 360, you will notice a circshift operation is necessary. This requires master branch versions of all packages I use, but I hope MakieOrg/GeoMakie.jl#29 will be resolved soon and we will get some stability.

here goes:

using Makie, GeoMakie, Proj4, MakieLayout
colorfunc(i) = circshift(A[:, :, i], 180)
Amean = circshift(drop(mean, A[:, :, 1:maxyearspan(t)], 3), 180) # average `A` over 19 years
texttime(i) = "$(attrib["long_name"]), units=[$units], t="*string(Date(t[i]))
titletext = "$(attrib["long_name"]), units=[$units], average of 19 years"

projection = "moll"
cmap = :inferno

source = Projection("+proj=lonlat +lon_0=180 +pm=180")
dest = Projection("+proj=$projection +lon_0=0")
xs, ys = xygrid(lons, lats)
Proj4.transform!(source, dest, vec(xs), vec(ys))

# get aspect ratio
xmin, xmax = extrema(xs)
ymin, ymax = extrema(ys)
aspect_ratio = (ymax - ymin) / (xmax - xmin)

total_crange = (minimum(A), maximum(A))
total_crange = (0, 400)

function plotfield(A, crange = (minimum(A), maximum(A));
                   cmap = :viridis, titletext)
    scene, layout = layoutscene(40; resolution = (1800, 950));
    # cmap = to_colormap(:viridis, 100)
    earthscene = layout[1, 1] = LScene(scene);

    sf = surface!(earthscene, xs, ys, zeros(size(xs)); color = A,
            shading = false, show_axis = false, colorrange = crange, colormap = cmap);

    geoaxis!(earthscene, -180, 180, -90, 90; crs = (src = source, dest = dest,));
    coastlines!(earthscene, 1; crs = (src = source, dest = dest,));

    colorbar = layout[1, 2] = LColorbar(scene, sf, width = 20);
    colsize!(layout, 1, Relative(1))
    rowsize!(layout, 1, Aspect(1, aspect_ratio))
    tt = layout[0, :] = LText(scene, titletext, textsize = 22);

    return scene, sf, tt
end

scene, sf, tt = plotfield(Amean; titletext = titletext, cmap = cmap)
Makie.save(plotsdir("ceresanim", "$(name)_avg.png"), scene)

# %% Then, make animation
scene, sf, tt = plotfield(A[:, :, 1], total_crange; titletext = titletext, cmap = cmap)
record(scene, plotsdir("ceresanim", "$(name).mp4"), 1:size(A, 3); framerate = 5) do i
    sf.color = colorfunc(i)
    tt.attributes.text[] = texttime(i)
end

this produces

and

Importantly, actually producing the animation requires as much time as looking at it, instead of the around 30 minutes it needs with cartopy.

[Balinus]

Really nice contribution (and figure)!

Do you want to do a PR? I can do it in the next couple of weeks on my end, probably around mid-march (some vacations I'm taking with the kids)

[gaelforget]

One possible caveat to adding geomakie as a dependency at this stage is that users may not be able to compile the package though -- at least on my Mac it fails miserably every time I try to use geomakie

see MakieOrg/GeoMakie.jl#29

[Datseris]

I can (and will) do a PR once we have a round of releases into all the packages that go into this process (they are like 6 of them!). At the moment @gaelforget has a very solid point. The good thing is that if there are working stable releases, we don't have to care about what breakages happen on master!

Although it is best to wait a bit because @asinghvi17 may create a simpler interface to using projections.

[asinghvi17]

@gaelforget could you post the error on GeoMakie? There will be some warnings, but I haven't seen any explicit failures to compile.

[gaelforget]

@gaelforget could you post the error on GeoMakie? There will be some warnings, but I haven't seen any explicit failures to compile.

Will do but it might take me a few days to get back to it. Might just have to do with some library that I don't know I need to install or something like that ...

[Balinus]

So, I made it work with one of the simulation I have on hand (CanESM2 CMIP5 simulation)! Yay!

I did get an error though on the coastlines!() function. Any ideas?

using following source and dest (my data is on a -180 to 180 reference.

source = Projection("+proj=lonlat +lon_0=0")
dest = Projection("+proj=$projection +lon_0=0")

coastlines!(earthscene, 1; crs = (src = source, dest = dest,));
ERROR: MethodError: no method matching coastlines!(::LScene, ::Int64; crs=(src = Projection("+proj=lonlat +lon_0=0 +ellps=GRS80"), dest = Projection("+proj=moll +lon_0=0 +ellps=GRS80")))
Closest candidates are:
  coastlines!(::Union{AbstractScene, AbstractPlotting.ScenePlot}; attributes...) at /home/proy/.julia/packages/GeoMakie/4QOPe/src/recipes/stock.jl:34
  coastlines!(::Attributes, ::Any...; kw_attributes...) at /home/proy/.julia/packages/GeoMakie/4QOPe/src/recipes/stock.jl:48
  coastlines!(::Union{AbstractScene, AbstractPlotting.ScenePlot}, ::Attributes, ::Any...; kw_attributes...) at /home/proy/.julia/packages/GeoMakie/4QOPe/src/recipes/stock.jl:53
Stacktrace:
 [1] #plotfield#19(::Symbol, ::String, ::typeof(plotfield), ::Array{Float32,2}, ::Tuple{Float32,Float32}) at ./REPL[196]:11
 [2] (::var"#kw##plotfield")(::NamedTuple{(:titletext,),Tuple{String}}, ::typeof(plotfield), ::Array{Float32,2}, ::Tuple{Float32,Float32}) at ./none:0
 [3] top-level scope at REPL[197]:1

[Balinus]

ok, it worked with the following call! (removed second positional argument)

coastlines!(earthscene; crs = (src = source, dest = dest,));

[Balinus]

Perhaps a simple question: how do we reverse the colormap?

[asinghvi17]

You can pass Reverse(:colormap) to the colormap kwarg :)

[Balinus]

Thanks! I knew there was something simple!

[Balinus]

I've read somewhere that Surface{...} is not supported by cairo. Is there any plan to add this feature? I am on a cluster without OpenGL capabilities. Thanks for any infos!

[asinghvi17]

Unfortunately there's no plan yet, but you can use meshes to simulate a surface. I can whip up an updated example later today.

[Balinus]

Yes, if I can fallback on meshes, that would be quite useful!

[gaelforget]

@gaelforget could you post the error on GeoMakie? There will be some warnings, but I haven't seen any explicit failures to compile.

Will do but it might take me a few days to get back to it. Might just have to do with some library that I don't know I need to install or something like that ...

Getting to this now -- sorry it's taken me longer than hoped to get things back on track after OSM20 ...

[asinghvi17]

Thanks, @gaelforget!

@Balinus, unfortunately Cairo meshes are not performant at all, and to generate even a 360x180 grid crashes CairoMakie. I'll look into this at Vizcon with the help of Simon, and perhaps work on a GR backend for Makie with the help of Josef Heinen, but till then if you don't have a GPU, you're out of luck. Will let you know if there's any more progress on that front.

[gaelforget]

@gaelforget could you post the error on GeoMakie? There will be some warnings, but I haven't seen any explicit failures to compile.

Will do but it might take me a few days to get back to it. Might just have to do with some library that I don't know I need to install or something like that ...

Getting to this now -- sorry it's taken me longer than hoped to get things back on track after OSM20 ...

Happy to report that I have been using GeoMakie v0.1.5 without any issue on MacOS & Julia v1.3.1 the past couple days. Not sure went wrong earlier...

One of my initial tests was this animation of ocean transports at 300m depth.

Next on my todo list would be to test GeoMakie on a binder instance -- has someone tried before?

[asinghvi17]

The problem with Binder is that it often doesn't have GPUs, which GeoMakie needs. You could try nextjournal.com? That offers free GPU compute.