Why would you want to achieve exactly the same thing twice but with a different language?
Although I use Plots.jl I personally don’t buy into the unified syntax story. Makie looks more promising in this regard because the backend abstraction is not based on existing languages but on hardware/drivers.
Maybe a rental analogy will do the trick: when you rent a car in the US or Europe, you don’t expect to have to pedal to open the trunk or to sing a local song to get the car started. It is reassuring that at least the basic controls are the same everywhere.
I agree. In its current iteration, the Julia “go-to” Plots package is amazing. 
What I do find difficult to follow is how all the plotting packages and backends fit with one another. In fact I’m so confused that I haven’t dared to go outside of Plots for the “publication-ready” stuff. (Oh and publication-ready is very much a subjective, field-dependent concept: some horrible plots do get published)
If someone in the know contributed a flow chart to illustrate which package to turn to according to need, it would greatly help our understanding of the Julia plotting universe. For instance, if you need a map, use X, Y. If you need 3D, use Y,Z. Etc.
There is such a table for Plots.jl backends. Backends · Plots
I’m sure most of those points apply for the packages the backends are based on too. The problem with comparing between frontends is most (advanced) users have committed. And even if you had a good objective comparison it’s likely to age poorly as the packages develop.
Personally I couldn’t speak to Makie, for instance, because I’ve dug down pretty hard in the Plots/GR camp.
@leon, It seems to be not the first time you are rising this topic. Is there any conclusion / summary of this and your previous conversations wrt your experience with Julia and Matlab? What do you think? IMO, even a few summary sentences would be interesting / useful in the near or long - term timeframe.
From my point of view, based on limited knowledge about all historical developments, I understand that Plots.jl was designed to somehow fulfill your suggestion. However, it seems that there was a will and a need for creation of other plotting packages that are not integrated with it. In general, I see it as a versatile world of Julia plotting. It has its advantages but also disadvantages. The main disadvantage for me is that I have to learn new things when I discover that some other package is providing abilities or ease of use that the one I am currently using is lacking. On the other side, I believe that such approach is moving things constantly forward. Having said it, I do believe that currently there is a way to solve a problem of “creating a figure with multiple subplots” when different packages are used, also at a level required for pro publications. At the same time, what I mentioned earlier, I would happily see a better support for basic geo plotting in the other packages than GMT and VegaLite. I believe for most of the people interested in geo plots, similar abilities as provided by VegaLite would be mostly sufficient. As addition, I will also allow myself to mention external package named kepler.gl that is not natively supported by Julia but is providing some abilities that are rather not present with any of Julia plotting packages currently. As a side note, I would like to strongly underline that I really had nothing wrong on my mind when addressing @carstenbauer a few posts above. Actually, I have to admit that I like his style of conversation very much, much better than sometimes suggested here so called niceties. 
Many thanks all for the great discussion.
I look forward to being able to plot figures like the below using GMT.jl in multiple subplots:
https://www.nature.com/articles/s41598-019-55039-4/figures/7
Also this one with one GMT.jl plotted map and another regular xy plot:
https://www.nature.com/articles/s41598-019-55039-4/figures/4
To start with, we must know the name of that projection. I though I knew but was wrong. Also, I’m not sure those plots were made with GMT. The use of the parula color scale make me think they were made with Matlab.
Search on the PROJ docs if you find that projection.
@Leon Just because those figures from the paper look nice I tried to replicate them with GeoMakie which is about to also have some major changes. For now, the following kinda works with scatters, (heatmaps and surfaces still are a problem, which I hope also get fixed in the next Makie release, don’t know really, @sdanisch probably will know better). Either way, here the code and the output 
using GeoMakie, CairoMakie
using Proj4, GeometryBasics, GeoDatasets
lonmask, latmask, data = GeoDatasets.landseamask(; resolution = 'c', grid = 10)
fig = Figure(resolution = (900, 1000), fontsize = 22)
ax1 = Axis(fig[1, 1], aspect = DataAspect(), title = "Projection: Goode H")
ax2 = Axis(fig[2, 1], aspect = DataAspect())
ax3 = Axis(fig[3, 1:2], xlabel = "date", ylabel = "y label")
# all this will go in the next release of GeoMakie
trans = Proj4.Transformation("+proj=longlat +datum=WGS84",
"+proj=igh_o +lon_0=-160", always_xy = true)
ptrans = Makie.PointTrans{2}(trans)
ax1.scene.transformation.transform_func[] = ptrans
ax2.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 = [1.15xytrans[i] for i in 1:prod(size(xytrans)) if checkInfs(xytrans[i])]
rectLimits = FRect2D(xytransClean)
limits!(ax1, rectLimits)
limits!(ax2, rectLimits)
# up to here things will go into the GeoMakie internals
field = [Point3f(lonmask[i], latmask[j], 0) for i = 1:3:length(lonmask) for j = 1:3:length(latmask)]
color1 = [Float64(data[i, j]) == 0.0 ? -rand() : Float64(data[i, j]) for i = 1:3:length(lonmask) for j = 1:3:length(latmask)]
color2 = [Float64(data[i, j]) == 0.0 ? -1 - rand() : Float64(data[i, j]) for i = 1:3:length(lonmask) for j = 1:3:length(latmask)]
mncolor = minimum(minimum.([color1, color2]))
mxcolor = maximum(maximum([color1, color2]))
hm = scatter!(ax1, field, color = color1, colorrange = (mncolor, mxcolor),
colormap = :CMRmap, markersize = 3)
scatter!(ax2, field, color = color2, colorrange = (mncolor, mxcolor),
colormap = :CMRmap, markersize = 3)
scatterlines!(ax3, 1 .. 10, x -> x^2)
Colorbar(fig[1:2, 2], hm, label = "some random label", width = 40, height = Relative(0.75))
hidedecorations!(ax1)
hidedecorations!(ax2)
hidespines!(ax1)
hidespines!(ax2)
save("GoodH.png", fig)
fig
Thank you so much for sharing! 
Over the years, I find the most effective way of learning a new tool is through a good working example. I’m sure a lot of people will find the sample code to be super useful.
Does GeoMakie allow a user to add GMT based plots to its subplots yet?
So I understand that GeoMakie is going to be in active development again, right?
Is the any estimated date for the next release?
Is GeoMakie going to be a separate project or is there a chance that GeoMakie will become a part of Makie with long term support?
I think it is called interrupted mollweide or something like that. It would be awesome if you could figure out how to plot such a map using GMT.jl. After all, GMT seems to be capable of producing the most vividly beautiful maps.
Below is my code to draw such a map in Matlab using the M_Map package. It plots the map piecemeal in 6 steps.
hold on;
set(gca,'xlimmode','auto','ylimmode','auto');
Slongs=[-100 43;-75 20; 20 145;43 100;145 295;100 295];
Slats= [ 0 90;-90 0;-90 0; 0 90;-90 0; 0 90];
for l=1:6,
m_proj('mollweide','long',Slongs(l,:),'lat',Slats(l,:));
m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...
'ytick',[-80:20:80],'yticklabels',[],'linest','-','color','k');
[C1, h1] = m_contourf (X, Y, Z, [5:0.01:10], 'LineStyle', 'none' );
[C1, h1] = m_contourf (X-360, Y, Z, [5:0.01:10], 'LineStyle', 'none' );
m_coast('patch',[.8 .8 .8]);
end;
set(gca,'xlimmode','auto','ylimmode','auto');
OK, so the problems. GMT has quite a lot o projections but nothing comparing to PROJ. It is in our plans to adopt PROJ for all the projections and in fact for point projections users can already use PROJ by using its syntax. Even for reprojecting images its is possible to use PROJ and all of its projections (not exactly true), but there is no mechanism to generate the frames and annotations for non-supported (by the GMT engine) projections. And the bad news is that the Interrupted Mollweide and Interrupted Good Homolosine are not supported projections. It should be possible to come out with a similar plot as shown by @lazarusA by plotting many points but I won’t go into the trouble of trying that.
A possible solution would be to try similar to what you do with M_Map and seeking inspiration in the Interrupted Sinusoidal example here
Ok, this is fun and interesting! I just took another look at the vega ecosystem, and I think with a few more tweaks to our Julia package we could make everything in https://github.com/d3/d3-geo-projection available as part of VegaLite.jl. And I think that also allows one to specify custom projections relatively easily. So my best guess is that one could take one of the projects in d3-geo-projection as a starting point (say geoInterruptedMollweide) and then customize it so that it looks like the example above.
If you are comfortable using LaTeX, I would honestly just make all your maps and subplots separately without axes. Then put them in a figure using pgfplots \addplot graphics. E.g. How do I fit a graphic exactly to an axis in pgfplots? - TeX - LaTeX Stack Exchange
That, or play around with saving subplots as images and doing something similar in Makie.
Well, we can still do nice plots like this but it’s not general enough to be very useful when data does not span on on land and oceans.
Here’s an example with bathymetry, but still struggling on how to manually make a frame.
using GMT
G = gmtread("@earth_relief_10m");
G2 = gdalwarp(G, ["-s_srs", "+proj=latlong", "-t_srs", "+proj=igh_o +lon_0=-160", "-dstnodata", "NaN"]);
imshow(G2, colorbar=true, shade=true, frame=:none)
It’s great to see GMT can do this type of plotting! Unfortunately, when I tested the code, I got something very different. See attached.
On a related question, does GMT allow subplots, so that I could plot 3 maps as shown in Figure 7 (https://www.nature.com/articles/s41598-019-55039-4/figures/7)? The projection does not have to be in interrupted Mollweide though.
You clearly plotted the geographic grid (G in my example) instead of the projected one (G2)
Yes. I even posted a link in one my posts in this thread.
It has also the very convenient options xshift and yshift
Does that use Jet as the default colormap? It has been replaced in many packages with “perceptually uniform” colormaps. See matplotlib colormaps for more info.
