Plot data on Map with GMT.jl

mixstam1453 · May 27, 2021, 3:11pm

Hi everyone,
I want to plot an 180x360 array on the Earth surface with coastlines. I know how to do it pycalling Cartopy, but i found GMT.jl interesting and usefull, but a little bit difficult to understand it as i am new on it.

1) plot data without coastlines , without correct lats!

A=rand(180,360)   # make the array
grdimage(A, fmt=:png,show=true) # plot it

2) plot data with coastlines, without correct lats!

grdimage(A)    
coast!(region=[0 360 -90 90],show=true)

3) plot only the HALF data in Robinson projection ( Northern Hemisphere) :

 A=rand(90,360) # grdimage doesn't want to plot A=rand(180,360)
topo = makecpt(color=:rainbow, range=(0.1,0.9,0.05), continuous=true);
grdimage(A,  coast=true,proj=:Robinson, frame=:a, color=topo,show=true)

######## 4) i want to plot my data in this figure:

coast(region=(-180, 180, -90, 90), 
    proj=(name=:Robinson, center=0),
    resolution=:crude,
    land=:white, water=:white, area=7000,figsize=15,
    shore=:thinnest,
    frame=(annot=45, grid=45,ticks=45),
    fmt = :jpg,show=true)

but i have no idea how to do it, in order to take all the latitudes correct…
Also, is there any way to do it with Plots.jl ? I didn’t find some way do it yet without errors, even with VegaLite, GeoMakie, GadFly, PlotlyJS…
Thanks for your time,
Michael

joa-quim · May 27, 2021, 6:06pm

Hi, thanks for trying it but won’t you prefer to plot something more colorful meaning than random noise? e.g.

grdimage("@earth_relief_01d", region=(-180,180,-90,90), coast=true, fmt=:png, show=true)

Note, you can’t make maps of Earth (or any other planetary body) without providing coordinates. rand(180,360) is just an array of numbers and grdimage just plotted in function of matrix size. To use your own data that does not come from known file formats (a very high number of them is recognized automatically) you must use the GMT.jl types. See docs for functions mat2grid, mat2img and mat2ds that help to quickly convert your data into the needed form.

See also Some examples · GMT for an example on how to plot data on maps. GMT itself has also many examples that look much nicer (and way less scaring) when using the GMT.jl wrapper.

mixstam1453 · May 28, 2021, 10:23am

Hi, congratulations for the package!
Thank you for your instructions.
I found that my problem was in mat2grid with incorrect x and y.
Finally i plot a 360x180 grided array of outgoing solar radiation, with Robinson projection:

B= Array(q[:,:,1]) # q is a ClimArray (by using ClimateBase), 180x360, but the lats are from 89.5 to -89.5
B=Array(B') # transpose the array in order to get 360x180
B=reverse(B, dims = 2)  # make the lats to be from -89.5 to 89.5
G2 = mat2grid(B , x=collect(0.5:1:359.5), y=collect(-89.5:1:89.5)) # x are the lons, y the lats

coast(region=(-180, 180, -90, 90), 
    proj=:Robinson,figsize=16,frame=(annot=45, grid=45,ticks=10)) # make the map
topo = makecpt(color=:rainbow, range=(0,225,20), continuous=true) #make colorbar
grdimage!(G2, coast=true,color=topo,
          frame=(annot=45,grid=45,title="CERES OSR")) # plot my gridarray
colorbar!(pos=(anchor=:BC,length=(12.5,0.6), horizontal=true, offset=(0,1.0)),
          color=topo, frame=(ylabel="W/m^2",), fmt=:jpg, show=true) # plot colorbar

Now i will looking for :
how to plot a text outside of the map (inside the map i do: text!([“Hello World”], x=92.0, y=30.0) ) ,
the color of title and text and
how to make the size of figure larger ( when figsize=30, some part of the right side of map is cutted…).
Thanks again!!

rafael.guerra · May 28, 2021, 11:43am

@mixstam, could you please edit your posts to include code inside ``` triple-backticks. See guidelines here. Thanks.

NB: are you using Linux or Windows?

mixstam1453 · May 28, 2021, 12:09pm

@ rafael.guerra Thank you. I use Linux, Kubuntu.

joa-quim · May 28, 2021, 12:57pm

OK, the image size issue first. This is something I must work more and at least document better. What happens is that by default we start with an A4 page size which the best for me that almost always do the visualization in PostScript. For all the other formats the page size is increased automatically but that it’s too late when the default output format is PS.
Short, for the moment, declare an environment variable called “JULIA_GMT_IMGFORMAT” and set it to "png". E.g., before starting GMT do

ENV["JULIA_GMT_IMGFORMAT"] = "png";

or after using GMT

GMT.FMT[1] = "png"

Now the text. Here the issue is that you want to plot the text somewhere in the universe. That is, outside the plot domain and that it’s not possible if you use the geographical coordinates. What you have to do is to start a new plot domain but this time using paper coordinates. Something like

text!(["Hello World"], region=(0,30,0,30), proj=:linear, x=15, y=15)

and you may have to add a figsize=... if the default one does not fit with previous plot.

mixstam1453 · May 31, 2021, 8:17am

Hello, have a good week.
Thanks again.
I plot my text well,

the code is here:


text(["mean = 170 W/m^2"],figsize=30,frame=(annot=""),region=(0,30,0,30), proj=:linear, x=7, y=9 ,font=(10,:red)) 
text!(["20 May 1994"],frame=(annot=""),region=(0,30,0,30), proj=:linear, x=1, y=9 ,font=(10,"Palatino-Italic",:darkgreen)) 
text!(["Outgoing Solar Radiation at TOA"],angle=7, frame=(annot=""),region=(0,30,0,30), proj=:linear, x=13, y=9.2 ,font=(10,"Courier-Bold",:purple))

coast!(region=(-180, 180, -90, 90), 
    proj=:Robinson,figsize=16,frame=(annot=45, grid=45,ticks=10)) # make the map
topo = makecpt(color=:rainbow, range=(0,225,20), continuous=true) #make colorbar
grdimage!(G2, coast=true,color=topo,
          frame=(annot=45,grid=45,title="CERES OSR")) # plot my gridarray
colorbar!(pos=(anchor=:BC,length=(12.5,0.6), horizontal=true, offset=(0,1.0)),
          color=topo, frame=(ylabel="W/m^2",), fmt=:jpg, show=true) # plot colorbar

but about the size of the image i can’t fix it with GMT.FMT[1] = "png" or ENV["JULIA_GMT_IMGFORMAT"] = "png"; , because when i execute this and then plot:

grdimage(G2, figsize=(30,14),coast=true,color=topo,frame=(annot=45,grid=45,title="CERES OSR")) # plot my gridarray

I get the larger size but when i try then colorbar!() or coast!() or text!() (something with !) i get this error:

julia> colorbar!(pos=(anchor=:BC,length=(12.5,0.6), horizontal=true, offset=(0,1.0)),             
                 color=topo, frame=(ylabel="W/m^2",), fmt=:jpg, show=true) # plot colorbar        
Error: /undefined in A
Operand stack:
   0
Execution stack:
   %interp_exit   .runexec2   --nostringval--   --nostringval--   --nostringval--   2   %stopped_push   --nostringval--   --nostringval--   --nostringval--   false   1   %stopped_push   1990   1   3   %oparray_pop   1989   1   3   %oparray_pop   1977   1   3   %oparray_pop   1833   1   3   %oparray_pop   --nostringval--   %errorexec_pop   .runexec2   --nostringval--   --nostringval--   --nostringval--   2   %stopped_push   --nostringval--
Dictionary stack:
   --dict:728/1123(ro)(G)--   --dict:0/20(G)--   --dict:75/200(L)--
Current allocation mode is local
Current file position is 9647080
psconvert [ERROR]: System call [gs -q -dSAFER -dNOPAUSE -dBATCH -sDEVICE=bbox -DPSL_no_pagefill -dMaxBitmap=2147483647 -dUseFastColor=true '/tmp/GMTjl_tmp.ps' 2> './psconvert_71425c.bb'] returned error 256.
ERROR: Something went wrong when calling the module. GMT error number =
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] gmt(::String)
   @ GMT ~/.julia/packages/GMT/wprG1/src/gmt_main.jl:289
 [3] showfig(d::Dict{Symbol, Any}, fname_ps::String, fname_ext::String, opt_T::String, K::Bool, fname::String)
   @ GMT ~/.julia/packages/GMT/wprG1/src/common_options.jl:2951
 [4] finish_PS_module(d::Dict{Symbol, Any}, cmd::Vector{String}, opt_extra::String, K::Bool, O::Bool, finish::Bool, args::GMT.GMTcpt)
   @ GMT ~/.julia/packages/GMT/wprG1/src/common_options.jl:3101
 [5] finish_PS_module(d::Dict{Symbol, Any}, cmd::String, opt_extra::String, K::Bool, O::Bool, finish::Bool, args::GMT.GMTcpt)
   @ GMT ~/.julia/packages/GMT/wprG1/src/common_options.jl:3028
 [6] colorbar(cmd0::String, arg1::Nothing; first::Bool, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{5, Symbol}, NamedTuple{(:pos, :color, :frame, :fmt, :show), Tuple{NamedTuple{(:anchor, :length, :horizontal, :offset), Tuple{Symbol, Tuple{Float64, Float64}, Bool, Tuple{Int64, Float64}}}, GMT.GMTcpt, NamedTuple{(:ylabel,), Tuple{String}}, Symbol, Bool}}})
   @ GMT ~/.julia/packages/GMT/wprG1/src/psscale.jl:85
 [7] #colorbar!#358
   @ ~/.julia/packages/GMT/wprG1/src/psscale.jl:93 [inlined]
 [8] top-level scope
   @ REPL[596]:1

Similar error i get even if i code:

text(["mean = 170 W/m^2"],figsize=30,frame=(annot=""),region=(0,30,0,30), proj=:linear, x=7, y=9 ,font=(10,:red)) 
text!(["20 May 1994"],frame=(annot=""),region=(0,30,0,30), proj=:linear, x=1, y=9 ,font=(10,"Palatino-Italic",:darkgreen))

julia> text!(["20 May 1994"],frame=(annot=""),region=(0,30,0,30), proj=:linear, x=1, y=9 ,font=(10,"Palatino-Italic",:darkgreen))                  
Error: /undefined in A
Operand stack:
   0
Execution stack:
   %interp_exit   .runexec2   --nostringval--   --nostringval--   --nostringval--   2   %stopped_push   --nostringval--   --nostringval--   --nostringval--   false   1   %stopped_push   1990   1   3   %oparray_pop   1989   1   3   %oparray_pop   1977   1   3   %oparray_pop   1833   1   3   %oparray_pop   --nostringval--   %errorexec_pop   .runexec2   --nostringval--   --nostringval--   --nostringval--   2   %stopped_push   --nostringval--
Dictionary stack:
   --dict:728/1123(ro)(G)--   --dict:0/20(G)--   --dict:75/200(L)--
Current allocation mode is local
Current file position is 20857
psconvert [ERROR]: System call [gs -q -dSAFER -dNOPAUSE -dBATCH -sDEVICE=bbox -DPSL_no_pagefill -dMaxBitmap=2147483647 -dUseFastColor=true '/tmp/GMTjl_tmp.ps' 2> './psconvert_71425c.bb'] returned error 256.
ERROR: Something went wrong when calling the module. GMT error number =
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] gmt(::String)
   @ GMT ~/.julia/packages/GMT/wprG1/src/gmt_main.jl:289
 [3] showfig(d::Dict{Symbol, Any}, fname_ps::String, fname_ext::String, opt_T::String, K::Bool, fname::String)
   @ GMT ~/.julia/packages/GMT/wprG1/src/common_options.jl:2951
 [4] finish_PS_module(::Dict{Symbol, Any}, ::Vector{String}, ::String, ::Bool, ::Bool, ::Bool, ::GMT.GMTdataset{Float64, 2}, ::Vararg{Any, N} where N)
   @ GMT ~/.julia/packages/GMT/wprG1/src/common_options.jl:3101
 [5] finish_PS_module(::Dict{Symbol, Any}, ::String, ::String, ::Bool, ::Bool, ::Bool, ::GMT.GMTdataset{Float64, 2}, ::Vararg{Any, N} where N)
   @ GMT ~/.julia/packages/GMT/wprG1/src/common_options.jl:3028
 [6] text(cmd0::String, arg1::GMT.GMTdataset{Float64, 2}; first::Bool, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:frame, :region, :proj, :font), Tuple{String, NTuple{4, Int64}, Symbol, Tuple{Int64, String, Symbol}}}})
   @ GMT ~/.julia/packages/GMT/wprG1/src/pstext.jl:127
 [7] text(txt::Vector{String}; x::Int64, y::Int64, first::Bool, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:frame, :region, :proj, :font), Tuple{String, NTuple{4, Int64}, Symbol, Tuple{Int64, String, Symbol}}}})
   @ GMT ~/.julia/packages/GMT/wprG1/src/pstext.jl:144
 [8] #text!#370
   @ ~/.julia/packages/GMT/wprG1/src/pstext.jl:146 [inlined]
 [9] top-level scope
   @ REPL[598]:1

To enlarge image size it is not my issue for now, but generally i would like to know how to achive this in future. I really appreciate your prompt responses and GMT.jl has a lot of features for further study.

Cheers,
Michael

joa-quim · May 31, 2021, 11:20am

Thanks. Opened an issue to track this. Diagnosis is done but need to come out with a solution.

joa-quim · May 31, 2021, 5:09pm

OK, the problem is hopefully solved but solution for now lives only on the master branch. It also now defaults to png so no more need to use the environment variable (unless you want to change the default to e.g. pdf)

Note, you can use region=:global and frame=:none instead of region=(-180, 180, -90, 90) and frame=(annot="")

mixstam1453 · June 1, 2021, 7:21am

Hi,that works !
Firstly ]add GMT#master
Here is the code with the larger figure size:

text(["mean = 170 W/m^2"],figsize=50,frame=:none,region=(0,50,0,50), proj=:linear, x=18, y=27 ,font=(18,:red)) 
text!(["20 May 1994"],frame=:none,region=(0,50,0,50), proj=:linear, x=37, y=27 ,font=(18,"Palatino-Italic",:darkgreen)) 
text!(["Outgoing Solar Radiation at TOA"],angle=-7, frame=:none,region=(0,50,0,50), proj=:linear, x=43, y=25.8 ,font=(18,"Courier-Bold",:purple))

coast!(region=:global, 
    proj=:Robinson,figsize=50) # make the map
topo = makecpt(color=:rainbow, range=(0,225,20), continuous=true) #make colorbar
grdimage!(G2,coast=true,color=topo,
          frame=(annot=30, grid=30,ticks=10,title="CERES OSR"),par=(FONT_ANNOT_PRIMARY=18,)) # plot my gridarray
colorbar!(pos=(anchor=:BC,length=(32.5,0.7), horizontal=true, offset=(0,1.0)),
          color=topo, frame=(ylabel="W/m^2",),par=(FONT_ANNOT_PRIMARY=18,), show=true,savefig="/home/michael/Desktop/gmt1.png") # plot colorbar & show and save figure

The figure size now is 50 and i also enlarge the font of annotations with par=(FONT_ANNOT_PRIMARY=18,)

Many thanks!
Cheers,
Michael

joa-quim · June 1, 2021, 11:28am

Cool and thanks for pushing me to do this fixes.

One thing on your plot. You are plotting the coastlines twice. First with the coast! command and then again with the coast=true option of grdimage. Furthermore the first coastlines plot is hidden by the next image plot. You may skip the call to coast! and move its options to grdimage where they are the same.

joa-quim · June 1, 2021, 3:28pm

@mixstam1453 A curiosity question. What is this Outgoing Solar Radiation at Top of Atmosphere? First I thought in ground reflection but looking at the image pattern it does not make much sense.

mixstam1453 · June 1, 2021, 3:53pm

Is the solar radiation that finally reflected back to the space (from clouds, aerosols and earth’s surface). The CERES project provides satellite-based observations of this. It is conceptually useful to think of the altitude at about 100 kilometers above the Earth as the “top of the atmosphere.”

joa-quim · June 1, 2021, 4:04pm

Thanks. On the back of my mind was a question on if this data could be useful to do better atmospheric correction of satellite images. But it would need to be on a per wavelength basis, not total radiation.

mixstam1453 · June 2, 2021, 7:41am

Yeah, i don’t study the satellite images but possibly it is necessary per wavelength basis for the correction.