Animation from grid data/files:

Specific Domains Visualization
,
1

Making an animation with the GMT.jl package is not intuitive to me and I would love some guidance/direction on how to achieve my end goal.

Luckily, though, using the GMT.jl package to plot a static grid is relatively straight forward, as it is only a few lines of code, shown below.

using DataFrames, DelimitedFiles, Glob, GMT,

epoch1 = readdlm("/home/rob/Documents/julia/Vgrdimg/51915_Vgrdimg.txt");

G = gridit(epoch1[:, [1, 2, 3]],
        mask = 0.3, V = :q);
grdimage(G, shade = true, coast = (ocean = :gray, shore = 0.5), fmt= :png, show, = true);

Which produced this image:

grid_gmt
grid_gmt1871×701 136 KB

And, it’s fine. Ideally, I would like the resolution of the gridding to be finer, like this plot I made in Python, but, I couldn’t get the masking to work in the eoMaps package.

sigma2_new
sigma2_new1000×500 222 KB

But I digress.

I tried to create a loop that would iterate through all my files, of which I have 1821, and each file is a timestep. So, I wanted to read the first 3-columns of the file, and then feed them into the movie() function in GMT, but I know I’m not doing it right, and I can’t figure out exactly what I’m supposed to do to get an animation to work.

Here is what I’ve attempted so far:

epoch1 = readdlm("/home/rob/Documents/julia/Vgrdimg/51915_Vgrdimg.txt");
file = "/home/rob/Documents/julia/Vgrdimg/"


function main_sc()
    for (i,file) in enumerate(glob("*.txt", file))
        file_hold = readdlm(file)

        G = gridit(file_hold[:, [1, 2, 3]],
        mask = 0.3, V = :q);

        fig = grdimage(G, shade = true, coast = (ocean = :gray, shore = 0.5), fmt= :png);

        #movie(main_sc, name=:count, frames = 100, format = :mp3, frame_rate = 16, clean = true );
        
    end
end
movie(main_sc, name=:count, frames = 100, format = :mp4, frame_rate = 16, clean = true );

Any help would be greatly appreciated. Thank you!

2

What do you get by calling

movie(main_sc, name=:count, frames = 100, format = :mp4, frame_rate = 16, clean = true );

and did you try gif?

3

Hi, I have not visited movie function for a while. What was the error?

Note, you can do this as well

epoch1 = gmtread("/home/rob/Documents/julia/Vgrdimg/51915_Vgrdimg.txt");

Then you must provide the inc=??? option to gridit (and probably adjust the value of the mask option). As is in your command the grid spacing wasw estimated from the data density.

This is not good (normally fatal) for movies. When you don’t provide the region it estimates it from data limits but there is no guarantie that the estimated limit is equal in all the frames. The solution is a fixed region=....

4 
epoch1 = gmtread("/home/rob/Documents/julia/Vgrdimg/51915_Vgrdimg.txt");

This is good to know! Thank you.

I did mess with various increments, but none of them actually made a, as far as I could tell, a difference in visualization. I tried 0.01, 0.1, 0.2, 0.5, 1.0, and 2.0 for the inc = option, and ultimately just left it out since it wasn’t making a difference.

Yeah, I didn’t think it would work, but I wanted to try something out at least. the GMT.jl documentations provides substantial detail for the input arguments, but, in the way of follow-able examples, it is certainly lacking.

I tried to understand and follow the example here because I think it’s a fairly close to what I’m trying to achieve.

When I run:

function main_sc()
    for (i,file) in enumerate(glob("*.txt", file))
        file_hold = readdlm(file)

        G = gridit(file_hold[:, [1, 2, 3]],
        mask = 0.3, V = :q);

        fig = grdimage(G, shade = true, coast = (ocean = :gray, shore = 0.5), fmt= :png);

        #movie(main_sc, name=:count, frames = 100, format = :mp3, frame_rate = 16, clean = true );
        
    end
end
movie(main_sc, name=:count, frames = 100, format = :mp4, frame_rate = 16, clean = true );

I get the following error message and stack trace:

ERROR: type String has no field data
Use `codeunits(str)` instead.
Stacktrace:
 [1] getproperty(x::String, f::Symbol)
   @ Base ./Base.jl:37
 [2] snif_GI_set_CTRLlimits(G_I::String)
   @ GMT ~/.julia/packages/GMT/PVGhC/src/imshow.jl:241
 [3] common_insert_R!(d::Dict{Symbol, Any}, O::Bool, cmd0::String, I_G::Nothing)
   @ GMT ~/.julia/packages/GMT/PVGhC/src/grdimage.jl:137
 [4] grdimage(cmd0::String, arg1::Nothing, arg2::Nothing, arg3::Nothing; first::Bool, kwargs::@Kwargs{…})
   @ GMT ~/.julia/packages/GMT/PVGhC/src/grdimage.jl:55
 [5] grdimage (repeats 2 times)
   @ ~/.julia/packages/GMT/PVGhC/src/grdimage.jl:51 [inlined]
 [6] main_sc()
   @ Main ~/Documents/julia/gridimage.jl:15
 [7] jl_sc_2_shell_sc(name::typeof(main_sc), name2::String)
   @ GMT ~/.julia/packages/GMT/PVGhC/src/movie.jl:164
 [8] movie(main::Function; pre::Nothing, post::Nothing, kwargs::@Kwargs{…})
   @ GMT ~/.julia/packages/GMT/PVGhC/src/movie.jl:110
 [9] top-level scope
   @ ~/Documents/julia/gridimage.jl:20

I also tried to see if I could decipher, modify, and translate this bash script to Julia syntax, but that is proving to be a job itself. lol

#!/bin/env -S bash -e

DIR0= $HOME/Documents/julia/Vgrdimg/
region=-R-130/-50/24.5/52
movie_files=movie_files
####### function to get the vaues of sigma max and min number of triangles from movie frame index
dir_names=($DIR0/*)
##~ get the name of the epoch from mnth_names that corresponds to the frame number
get_values() {
    local index=$1
    echo "${mnth_names[index]}"
}

out=strain_map_year ##not sure this is even needed 
gmt begin "$out"
    ###########################~ plot coastlines, strain grid   ######################
    gmt coast $region -B --FORMAT_GEO_MAP=dddF -JM15c -EUS.CA,US.OR+p0.5p,#6e04c4 -W0.1p,#b4b5cc -G#ddded9 -S#b4b5cc -Df
    gmt grdimage "$movie_files"/ave_strain_${MOVIE_FRAME}.grd -Cdil_strain.cpt
    ###########################~ plot principal strains  #############################
    pstr_ampl=0.08
    ##############& e1 is compressional
    awk -v f="$pstr_ampl" '{if($3<0) print $1, $2, $5, $3*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.75c+ea+je+n1/0.2 -W3p,black -Gblack #-SVh0.04c/0.085c/0.075c -W0.1,black -Gblack
    awk -v f="$pstr_ampl" '{if($3<0) print $1, $2, $5-180, $3*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.75c+ea+je+n1/0.2 -W3p,black -Gblack
    ###############& e2 is compressional
    awk -v f="$pstr_ampl" '{if($4<0) print $1, $2, $5+90, $4*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.75c+ea50+je+n1/0.2 -W3p,black -Gblack
    awk -v f="$pstr_ampl" '{if($4<0) print $1, $2, $5-90, $4*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.75c+ea+je+n1/0.2 -W3p,black -Gblack
    ################& e1 is dilatational
    ######* the first 2 lines are there to have a black outline:
    awk -v f="$pstr_ampl" '{if($3>0) print $1, $2, $5, $3*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.75c+n1/0.2 -W2p,black -Gblack
    awk -v f="$pstr_ampl" '{if($3>0) print $1, $2, $5-180, $3*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.75c+n1/0.2 -W2p,black -Gblack
    ###### *the second 2 lines actually show the sizes of the strain axes in white
    awk -v f="$pstr_ampl" '{if($3>0) print $1, $2, $5, $3*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.7c+ea+jb+n1/0.2+p0.4p -W0.5p,white -Gwhite
    awk -v f="$pstr_ampl" '{if($3>0) print $1, $2, $5-180, $3*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.7c+ea+jb+n1/0.2+p0.4p -W0.5p,white -Gwhite

    ###################& e2 is dilatational
    ######* the first 2 lines are there to have a black outline:
    awk -v f="$pstr_ampl" '{if($4>0) print $1, $2, $5+90, $4*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.75c+n1/0.2 -W2p,black -Gblack
    awk -v f="$pstr_ampl" '{if($4>0) print $1, $2, $5-90, $4*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.75c+n1/0.2 -W2p,black -Gblack
    ###### *the second 2 lines actually show the sizes of the strain axes in white    
    awk -v f="$pstr_ampl" '{if($4>0) print $1, $2, $5+90, $4*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.7c+ea+jb+n1/0.2+p0.4p -W0.5p,white -Gwhite
    awk -v f="$pstr_ampl" '{if($4>0) print $1, $2, $5-90, $4*f}' $DIR0/MELD/seaspred/"$movie_files"/principal_strains${MOVIE_FRAME}.gmt | gmt plot -SV0.7c+ea+jb+n1/0.2+p0.4p -W0.5p,white -Gwhite #-SVt0.04c/0.085c/0.075c -W0.1,black -Gwhite

########################~ plot the name of epoch, sigma_max, and Nmin  ##############
    month_vals=$(get_values ${MOVIE_FRAME})
    gmt text -F+f17p,Helvetica-Bold+jCB << EOF
        -119 45.8 sigma: 25, Nmin: 40
        -122 45.8 $month_vals
EOF
########################~ add stateline, lake outlines on top of frid  ##############
    gmt coast -W1 -Df -Na -Lf-130.8/46/10/200+lkm -EUS.CA+p0.5p
###########################~   plot mountain(s)   ###################################
    gmt plot mthood.lonlat -St0.5c -Wyellow -Gblack
###########################~  plot med locs of seism  ##########################
    gmt plot extra_locations.txt -Sc0.5c -Wblack -G#e00757
###########################~ plot seismic swarms   ##################################
    # sed 's/,/\t/g' maupin_swarm.csv | awk '{print $1, $2, $3, $4}' > maupin_swarm_scol.txt #date, lat, lon, depth, mag 
    ## by depth
    # awk '{print $3, $2, $4, $5}' $DIR0/seismo_data/maupin_swarm_scol.txt | gmt plot -Cseismo_depth.cpt -Sci -i0,1,2,3+s0.06 -Wfaint -hi1
    # gmt psscale -Cseismo_depth.cpt -DjBR+o-1.5c/-1.2c+w6c/0.25c+h+o0c/0c -Baf+l"Seismic depth (km)" -I --FONT_LABEL=14p
    ## by date
    ###& Maupin swarm
    awk '{print $5, $4}' $DIR0/seismo_data/maupin_swarm_locations_decyr_IEB.txt | gmt plot -Sc0.25c -hi1 -G#f7bc0a -Wblack
    ###& NLT swarm 
    awk '{print $3, $2, $1}' $DIR0/seismo_data/north_lake_tahoe_2003_swarm_decyr.txt | gmt plot -Sc0.25c -hi1 -G#f7bc0a -Wblack
    ###& Sierra Valley swarm
    awk '{print $5, $4}' $DIR0/seismo_data/sierra_valley_locations_decyr_trugman.txt | gmt plot -Sc0.25c -hi1 -G#f7bc0a -Wblack
    # gmt psscale -DjBR+o-1.5c/-1.2c+w6c/0.25c+h+o0c/0c -Baf+l"EQ year" -I --FONT_LABEL=14p -Cseismo_date.cpt 
    ###& focal mech(s)
    gmt psmeca $DIR0/seismo_data/single_maupin_moment_tensor.gmt -Sa1c -A

###########################~ add scale for strain  ###################################
    gmt psscale -Cdil_strain.cpt -DjBL+o4c/-1.2c+w6c/0.25c+h+o0c/0c -Baf+l"Dilatational Strain (10^-9)" -I --FONT_LABEL=14p
gmt end show


#### use these the make movie from command line:

####! use these the make movie from command line:
##& make a single frame to look at (here i have 24 frames, so -T24; -M is the master frame to view to make sure it looks right (as a png), so -M6 makes the 6th frame; -Fnone tells it we're not making a full movie yet); -Z deletes the intermediate directory with the files needed to make the movie; -N is the name of the output frame or output movie; -C sets the dimensions of the canvas; -D is how long you want each frame to display (frame rate))
##~ make 6th frame
#  gmt movie make_strain_movie_fullyear.sh -C20x25x200 -T24 -M6,png -Fnone -Nseason_strain -Z

##~ make actual movie: 
##&-P sets the options for a progress indicator (I have a bar on top right with a red triangle that shows where in the movie you are)
# gmt movie make_strain_movie_fullyear.sh -C20x25x200 -T24 -Fmp4 -Nseason_strain -Z -D0.5 -Pf+jTR+o0.5
5

I saw the option for gif but, I wasn’t going to mess around with option arguments until I get the required ones implemented correctly.

Running this:

using DataFrames, DelimitedFiles, Glob, GMT

epoch1 = readdlm("/home/rob/Documents/julia/Vgrdimg/51915_Vgrdimg.txt");
file = "/home/rob/Documents/julia/Vgrdimg/"


function main_sc()
    for (i,file) in enumerate(glob("*.txt", file))
        file_hold = readdlm(file)

        G = gridit(file_hold[:, [1, 2, 3]],
        mask = 0.3, V = :q);

        fig = grdimage(G, shade = true, coast = (ocean = :gray, shore = 0.5), fmt = :png);

        #movie(main_sc, name=:count, frames = 100, format = :mp3, frame_rate = 16, clean = true );
        
    end
end
movie(main_sc, name=:count, frames = 100, format = :mp4, frame_rate = 16, clean = true );

I get this error message:

ERROR: type String has no field data
Use `codeunits(str)` instead.
Stacktrace:
 [1] getproperty(x::String, f::Symbol)
   @ Base ./Base.jl:37
 [2] snif_GI_set_CTRLlimits(G_I::String)
   @ GMT ~/.julia/packages/GMT/PVGhC/src/imshow.jl:241
 [3] common_insert_R!(d::Dict{Symbol, Any}, O::Bool, cmd0::String, I_G::Nothing)
   @ GMT ~/.julia/packages/GMT/PVGhC/src/grdimage.jl:137
 [4] grdimage(cmd0::String, arg1::Nothing, arg2::Nothing, arg3::Nothing; first::Bool, kwargs::@Kwargs{…})
   @ GMT ~/.julia/packages/GMT/PVGhC/src/grdimage.jl:55
 [5] grdimage (repeats 2 times)
   @ ~/.julia/packages/GMT/PVGhC/src/grdimage.jl:51 [inlined]
 [6] main_sc()
   @ Main ~/Documents/julia/gridimage.jl:14
 [7] jl_sc_2_shell_sc(name::typeof(main_sc), name2::String)
   @ GMT ~/.julia/packages/GMT/PVGhC/src/movie.jl:164
 [8] movie(main::Function; pre::Nothing, post::Nothing, kwargs::@Kwargs{…})
   @ GMT ~/.julia/packages/GMT/PVGhC/src/movie.jl:110
 [9] top-level scope
   @ ~/Documents/julia/gridimage.jl:1
Some type information was truncated. Use `show(err)` to see complete types.
6

Strange, it should. Can you make available one of those 51915_Vgrdimg.txt for me too if something is nor working as it should.

I’ll be looking at the movie issues. Maybe it’s better to continue this in the GMT forum where another user is also asking movie related questions.

7

Sure thing. Oh, guess I can’t attach text files. Let I’ll try to upload the files to my cloud storage and share that link.

Here’s the link

The file values are as follow:

[lon, lat, value, std_dev, std_dev_alt, num_stations_used, ave_distance(km)]

And I’m only plotting the first 3 columns.

I’ll do that, thanks for the suggestion and for your help with my problem!

8

inc = is working now, and the resolution has greatly improved. I’m not sure what was going on yesterday. Mostly likely I had it somewhere in the operation where It wasn’t getting read.

GMT_roma
GMT_roma1972×1119 367 KB

 # value for cpt in epoch1[:, 3]
C = GMT.makecpt(cmap=:roma, 
    range=(minimum(epoch1[:, 3]), maximum(epoch1[:, 3]), 0.075));
#imshow(C, B = :none, horizontal = false) # check increments of colorbar

G = GMT.gridit(epoch1[:, [1, 2, 3]], mask=0.3, inc = 0.1, V=:q);


GMT.grdimage(G, 
    proj=(name = :lambertConic, center = [-100 35], parallels = [33 45]),
    coast = (ocean = :gray, shore = 0.5),
    frame=:agf,
    cmap=C, 
    title = "Vertical Residual Values", 
    colorbar = true,
    fmt = :png, 
    show = true 
    );

But this is much better. Now, just to figure out how to animate it…

9

Good that it worked (one less thing to look at). As I said I’ll try to see the animation part. Some capabilities of the GMT CLI movie are hard do reproduce in GMT.jl because the C lib generates and executes shell scripts for intermediate steps but I think this example should work.

On a side note, you can simplify a little your commands by doing:

epoch1 = gmtread("51915_Vgrdimg.txt", incols="0,1,2")
BoundingBox: [-124.4, -53.0, 24.8, 49.2, -6.794, 4.959]

31104×3 GMTdataset{Float64, 2}
   Row │  col.1  col.2   col.3
───────┼───────────────────────
     1 │  -82.6   24.8  2.8215
     2 │  -82.4   24.8  2.8215
     3 │  -82.2   24.8  2.8215
     4 │  -82.0   24.8  2.8215

G = gridit(epoch1,  mask=0.3, inc = 0.1);
2 Likes
10

Can you call the following each time you load a txt file:

epoch1 = gmtread("51915_Vgrdimg.txt", incols="0,1,2")
@assert typeof(epoch1) == GMTdataset{Float64, 2}
@show epoch1
@show typeof(epoch1)
11

Hi, I have opened this GMT.jl issue that explains some of the details involved in the movie usage and to track future improvements.

But regarding your case, I think that if you rename the input .txt files to *_Vgrdimg_N.txt, where N = 0:n_files-1, the next piece of code should a working base that you can tune. Note, you must adapt the frames = ??? to your number of files.

function pre_sc()
	makecpt(cmap=:roma, range=(-6.8, 5.0, 0.1), H=true, name="lixo.cpt");
end
function main2_sc()
	surface("51915_Vgrdimg_" * "MOVIE_FRAME" * ".txt", incols="0,1,2", proj=(name=:lambertConic, center=[-90 35], parallels=[33 45]), region=(-124.4, -53.0, 24.8, 49.2), inc=0.1, mask=0.3, save="lixo.grd");
	grdimage("lixo.grd", shade=true, cmap="lixo.cpt", proj=(name=:lambertConic, center=[-90 35], parallels=[33 45]), frame=:auto, title = "Vertical Residual Values")
	coast!(ocean=:gray, shore=0.5)
end
movie(main2_sc, pre=pre_sc, C="hd", name=:deforms, frames=10, format=:mp4, frame_rate=10, V=true)

Put this in a script and include(it.jl)

EDIT: Sorry, forgot to say. Before a version > 1.13.3 comes out, the master version is needed to run the above.

12

If the movie is the last sticking point, then a workaround could be (assuming your’e on linux based on the filepaths) to use ffmpeg directly. For a series of images, such as image_001.png, image_002.png, …, it is easy to make a movie by doing (assuming no typos)

ffmpeg -framerate 16 -i 'image_*.png' -o movie.mp4

If desired, it might then be possible to pop that into your Julia script and do

run(`ffmpeg -framerate 16 -i 'image_*.png' -o movie.mp4`)

Just an option!

1 Like
13

Yes, that’s true. And we even mention such option in the GMT CLI man page.

1 Like