Animation : sea-ice concentration in the Arctic Ocean

Hi there !

Module used :

  • mapproject , grdproject , math , grdsample
  • makecpt , colorbar
  • grdimage , coast (including clipping) , basemap
  • clip
  • movie

Full script :

Animation : sea-ice cover in Arctic sea (2016).
Stereographic projection (-Js) with some regions leaking out of the frame
(based on IBCAO map)
DATASET :
topography and bathymetry
sea-ice

Include.sh

lon_ref=0           # map centered on Greenwich meridian
lat_ref=90          # map centered on North pole

lat_min=65          # map extend to 65th parallele
lat_out=55          # (extra) map extend to cover Greenland
lat_scl=77          # true scale at this parallele
scale=1:40000000    # map size (here ~ 14cm)

# Calculation used to properly place the plotted regions falling out of the frame
W1=$(gmt mapproject -R0/360/${lat_min}/90 -Js${lon_ref}/${lat_ref}/${lat_scl}/${scale} -Ww)
H1=$(gmt mapproject -R0/360/${lat_min}/90 -Js${lon_ref}/${lat_ref}/${lat_scl}/${scale} -Wh)
W2=$(gmt mapproject -R0/360/${lat_out}/90 -Js${lon_ref}/${lat_ref}/${lat_scl}/${scale} -Ww)
H2=$(gmt mapproject -R0/360/${lat_out}/90 -Js${lon_ref}/${lat_ref}/${lat_scl}/${scale} -Wh)
	
dx=$(gmt math -Q $W1 $W2 SUB 2 DIV =)
dy=$(gmt math -Q $H1 $H2 SUB 2 DIV =)

Pre.sh


gmt begin
    gmt makecpt -Cbukavu -T-4500/1500 -H > relief.cpt                                                   # Colormap for topography
    gmt makecpt -Cvik -G-0.5/0 -T0/101/1 -H -N > ice.cpt                                                # Colormap for sea-ice

    gmt math -o0 -fT -T2016-01-01/2016-12-31/1 T --TIME_UNIT=d --FORMAT_DATE_OUT=yyyymmdd = list.txt    # List of file's extension used as movie-frame (yyyymmddTHH:MM:SS)
    awk -F[T] '{print $1}' list.txt > list2.txt                                                         # Trim to yyymmdd

    gmt grdimage -R0/360/${lat_min}/90 -Js${lon_ref}/${lat_ref}/${lat_scl}/${scale} @earth_relief_01m -Crelief.cpt -I+d -Y4c    # DEM inner region with illumination, shifted up for colorbar
    gmt coast -Dh -Wthinnest,black -Ir/faint,lightblue  -Baf                                                                    # corresponding coast lines (+rivers) with framebox

    gmt coast -ECA.NU,GL,IS+c -R0/360/${lat_out}/90 -X${dx} -Y${dy}                                                             # Prepare outter region for clipping (shifted to match first layers)
    gmt clip -N << END                                                                                                          # ... but exclude this box
-80.75 ${lat_min} 
-120 ${lat_min} 
-120 ${lat_out}
-80.75 ${lat_out}
END
    gmt grdimage @earth_relief_01m -Crelief.cpt -I+d                                                                            # DEM outter region with illumination
    gmt coast -Dh -Wthinnest,black -Ir/faint,lightblue                                                                          # corresponding coast lines (+rivers)
    gmt clip -C                                                                                                                 # close excluded box clipping
    gmt coast -Q                                                                                                                # close outter region clipping

    gmt colorbar -Cice.cpt -DJMR+w${H1}c/0.5c+o-1.5c/0c -G0/100 -Np -Bxa10f1g10+l"Sea-ice concentration" -By+l"%"               # Colorbar (right) for sea-ice
    gmt colorbar -Crelief.cpt -DJBC+w${W1}c/0.5c+o0c/-1.5c+h+e -G-4000/1000 -Bafg+l"Topography" -By+l"m" -I                     # Colorbar (bottom) for relief
gmt end

Main.sh

gmt begin

    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
    # Operation on SEA-ICE coverage data :                                          #
    #   1. Reproject from (x,y) [km relative to the pole] to (lon,lat) [degree]     #
    #   2. Upsample using bilinear interpolation AND remove zeros (=NaNs)           #
    # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

    gmt grdproject -R0/360/${lat_min}/90 -Js${lon_ref}/${lat_ref}/${lat_scl}/${scale} data/data_${MOVIE_COL0}.nc?"ice_conc" -Gdata_bis.nc -C -I -Fk
    gmt grdsample data_bis.nc+n0 -nl -I1m -Gdata_ter.nc

    gmt grdimage -R0/360/${lat_min}/90 -Js${lon_ref}/${lat_ref}/${lat_scl}/${scale} data_ter.nc -Cice.cpt -Q  -Y4c                              # plot sea-ice (ignore NaN's), aligned with main frame
    gmt basemap -B+t"$(gmt math -Q -o0 -fT ${MOVIE_COL0}T --FORMAT_DATE_IN=yyyymmdd --FORMAT_DATE_OUT=yyyy-o-dd = | awk -F[-] '{print $2,$1}')" # title ([abbr. month]<space>[year])

gmt end

# Movie command (with progress circle)
gmt movie main.sh -Iinclude.sh -Sbpre.sh -Narctic_mov -C20cx21.5cx120 -Tlist2.txt -D18 -Fmp4 -Pb+jTL+o1c/1.25c

(the glitch in March is due to a missing file. I should have check more carefully)

1 Like

Impressive!! Very nice!!

1 Like

Thanks :slight_smile:
(But so far nothing can compete with your Messi animation)

1 Like

Really slick. Would you be interested in making an animation for inclusion to the animation gallery? I.e., make it HD and ideally pull the sea-icea data from a data place (we would like not having to store the data grids)?

It is HD no ?
But sure, I can make that… Though, the data can be heavy and/or long to download. What do you suggest ?

edit : FYI, 1 year of this dataset (daily, gridded) is ~3.6GB

By “HD” I mean 1980x1280 in 16:9 format. I think it is important that all our animations are reproducible, even if they take a long time to run or download. For instance, some of the animations takes a few hours to render, so I would not worry about a big download for those who wish to try and modify. I am just hoping there are not too many custom processing steps you did. Ideally, it is just a data download as part of pre.sh and then off to build the movie.

So you may need to make a few adjustments to fit the 24 x 13.5 cm HD “papersize”

Maybe two of these side by side, or Arctic vs Antarctic?

My nemesis ! The -C option of movie xD

The only thing I did to the dataset was to rename it to a shorter format (data_20160101.nc instead of ESACCI-SEAICE-L4-SICONC-AMSR_25.0kmEASE2-SH-20160101-fv2.1.nc)

Everything else is “as is” in the first post.

For the download from ECMWF, an account is required (a token is given for the API, then you can use a python script to download everything)

The subplot module doesn’t like -X, -Y shift … so it’s going to be just one big figure with a big shift between the two plots (-X15c or something)

Yes, may need to be a single figure with one part shifted by -X.

[updated video next post]

There was a lot of errors along the process. GMT is not yet very stable when the load is distributed on many cpus (24 here).

Here are some examples of the output I had (corresponding to the video uploaded)

bash poles_movie.sh 

movie_init.sh: line 29: ./gmt_project.sh: Text file busy
grdproject [ERROR]: Cannot find file 3
grdsample [ERROR]: Cannot find file data_bis.nc
gmtmath [ERROR]: Operation "ADD" requires 2 operands
grdproject [WARNING]: For a UTM or TM projection, your region -1.000000000000/1.000000000000/89.000000000000/91.000000000000 is too large to be in degrees and thus assumed to be in meters
grdsample [WARNING]: (y_max-y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.
grdsample (gmtapi_init_grdheader): Please select compatible -R and -I values
grdsample [WARNING]: Output sampling interval in y exceeds input interval and may lead to aliasing.
grdimage [ERROR]: Must specify one (or three [deprecated]) input file(s)
gmtmath [ERROR]: Operation "ADD" requires 2 operands
[...]

Also, for some reason the basemap for Antarctica doesn’t display the annotation…

[update : I’m downloading 2013, 2014 and 2015 to have a longer animation with a higher frame rate]

1 Like

@pwessel

Longer version (2013 - 2016) :

Script :

(someone more seasoned should review it for improvement, especially :

  • the “cp + chmod” in include.sh)
  • the list creation in pre.sh
  • the file management (download - unzip - sort)
  • add the -W option to work in a more adapted directory
  • maybe add the screen command (I use it but it’s not in the script)?
  • I was also thinking about adding a graph at the bottom showing the sea-ice coverage as a function of time … If someone is motivated enough to implement that :slight_smile:

poles_movie.sh.txt (11.9 KB) [remove .txt extension]]

(nota : the video is nicer at 0.5x speed…)

1 Like

Thanks @PlanetGus I will have a look as soon as I can. I am also having some issues with 24 cores when movie launches jogs that deal with curl…

Regarding the “ice-coverage” plot I mentioned adding …
Do you have a more clever use of grdmath ? Creating a 2-D grid for 1 value doesn’t seem efficient.

projection_north='-Js0/90/77/1:93300000'
region_north='-R0/360/50/90'

for ff in $(ls -v data/*nc);
do 

dd=$(ls ${ff} | awk -F[_.] '{print $2}')
gmt grdproject $projection_north $region_north ${ff}?ice_conc -Gdata_bis.nc -C -I -Fk
gmt grdmath data_bis.nc 1 GE SUM  = tmp.nc
gmt math -Q tmp.nc?z[0] = | awk -v day=${dd} 'NR==1 {print day,$0}' >> x.txt
done

gmt begin ice_cov png 
    gmt plot -R2016-01-01T/2016-12-31T/0/150000 -JX10cT/10c x.txt -Wred --FORMAT_DATE_IN=yyyymmdd -B+t"Sea-ice coverage" -l"North Pole" -Bxafg -Byafg+l"cell count"
gmt end show

(Also, converting cell count to km^2 could be a nice thing to do)

Sorry for not following up - will soon - but swamped with a few things.