Example scripts for eddy composites over Agulhas rings

eddy_composite/IDeddy_mm.py

@@ -0,0 +1,61 @@
+#!/bin/python
+
+# ©2023 MPI-M, Dian Putrasahan
+
+import sys
+yyyy=sys.argv[1]
+mth=sys.argv[2]
+
+print('Identifying monthly eddies for year='+str(yyyy)+', mth='+str(mth))
+
+from py_eddy_tracker.dataset.grid import RegularGridDataset
+from datetime import datetime
+
+#Read in example SSH data that has been mapped onto a 0.25deg regular grid.
+# /work/mh0287/m300466/topaz/ngc2013/zos
+# /work/mh0287/m300466/topaz/rthk001/zos
+
+expid='ngc2013'
+varname='zos'
+fq='mm'
+outdir='/work/mh0287/m300466/topaz/'+expid+'/'+varname+'/'+fq+'/'
+date = datetime(int(yyyy), int(mth), 1)
+grid_name, lon_name, lat_name = (
+    outdir+expid+'_'+varname+'_'+fq+'_'+date.strftime('%Y%m')+'_MR25.nc',
+    "lon",
+    "lat",
+)
+g = RegularGridDataset(grid_name, lon_name, lat_name) #assumes no time index, might only take the first time index
+
+from py_eddy_tracker import data
+from py_eddy_tracker.dataset.grid import RegularGridDataset
+
+# from py_eddy_tracker import start_logger
+# start_logger().setLevel("DEBUG")  # Available options: ERROR, WARNING, INFO, DEBUG
+
+##Load Input grid, SSH is used to detect eddies. Add a new filed to store the high-pass filtered SSHA
+g.add_uv("zos")
+g.copy("zos", "zos_high")
+wavelength = 700  #choice of spatial cutoff for high pass filter in km
+g.bessel_high_filter("zos", wavelength, order=1)
+
+# Run the detection for the total grid and the filtered grid
+step_ht=0.005 #intervals to search for closed contours (5mm in this case)
+a, c = g.eddy_identification(
+    "zos_high", "u", "v",
+    date,  # Date of identification
+    step_ht,  # step between two isolines of detection (m)
+    pixel_limit=(50, 400),  # Min and max pixel count for valid contour
+    shape_error=30,  # Error max (%) between ratio of circle fit and contour
+)
+
+#Save output
+from netCDF4 import Dataset
+
+with Dataset(date.strftime(outdir+expid+"_anticyclonic_"+fq+"_"+date.strftime('%Y%m')+".nc"), "w") as h:
+    a.to_netcdf(h)
+with Dataset(date.strftime(outdir+expid+"_cyclonic_"+fq+"_"+date.strftime('%Y%m')+".nc"), "w") as h:
+    c.to_netcdf(h)
+
+
+

eddy_composite/LICENSE

@@ -0,0 +1,28 @@
+BSD 3-Clause License
+
+Copyright (c) 2023, EERIE 
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

eddy_composite/README.md

@@ -0,0 +1,46 @@
+# eddy_composites
+
+Scripts and notebooks for building eddy composites with example over Agulhas rings
+©2023 MPI-M, Dian Putrasahan
+
+1) Remap data onto MR25 grid
+	cd /home/m/m300466/NGC2/eddytrack
+	>> sbatch remap_mm_r2b8O_MR25.job zos ngc2013
+
+2) Split files into yyyymm 
+	Edit  split_mm_yrmth.sh
+	varname=zos
+	varname=Wind_Speed_10m
+	varname=atmos_fluxes_HeatFlux_Latent
+	varname=mlotst
+	varname=sea_level_pressure
+	varname=to
+	expid=ngc2013
+	expid=rthk001
+	>> ./split_mm_yrmth.sh
+
+3) Eddy identification to get anticyclones and cyclones (only with zos)
+	cd /home/m/m300466/NGC2/eddytrack
+	Notebook with some figures to show identified eddies [eddy-identification.ipynb]
+	#Submit scripts for each year using submit_IDeddy_mm.job, which in turn submits IDeddy_mm.py script for each month
+	>> for yr in $(seq 2020 2049); do sbatch submit_IDeddy_mm.job ${yr}; done
+
+4) Remap and high pass filter other fields
+	cd /home/m/m300466/NGC2/eddytrack
+	>> sbatch remap_mm_r2b8O_MR25.job Wind_Speed_10m ngc2013
+	>> sbatch remap_mm_r2b8O_MR25.job atmos_fluxes_HeatFlux_Latent ngc2013
+	>> sbatch remap_mm_r2b8O_MR25.job mlotst ngc2013
+	>> sbatch remap_mm_r2b8O_MR25.job sea_level_pressure ngc2013
+	>> sbatch remap_mm_r2b8O_MR25.job to ngc2013
+	# Submit scripts for each year using submit_highpass_mm.job, which in turn submits highpass_mm.py script
+	>> for yr in $(seq 2020 2049); do sbatch submit_highpass_mm.job zos ${yr}; done
+	>> for yr in $(seq 2020 2050); do sbatch submit_highpass_mm.job Wind_Speed_10m ${yr}; done
+	>> for yr in $(seq 2020 2050); do sbatch submit_highpass_mm.job atmos_fluxes_HeatFlux_Latent ${yr}; done
+	>> for yr in $(seq 2020 2050); do sbatch submit_highpass_mm.job to ${yr}; done
+	>> for yr in $(seq 2020 2050); do sbatch submit_highpass_mm.job atmos_fluxes_FrshFlux_Precipitation ${yr}; done
+	>> for yr in $(seq 2020 2050); do sbatch submit_highpass_mm.job mlotst ${yr}; done
+
+5) Composite for a given region (Agulhas rings)
+monthly-eddy-composites.ipynb
+
+

eddy_composite/highpass_mm.py

@@ -0,0 +1,39 @@
+#!/bin/python
+# ©2023 MPI-M, Dian Putrasahan
+
+import sys
+varname=sys.argv[1]
+yyyy=sys.argv[2]
+mth=sys.argv[3]
+#varname='zos'
+#varname='Wind_Speed_10m'
+
+print('High pass filter monthly '+varname+' for year='+str(yyyy)+', mth='+str(mth))
+
+from py_eddy_tracker.dataset.grid import RegularGridDataset
+from datetime import datetime
+
+expid='ngc2013'
+fq='mm'
+outdir='/work/mh0287/m300466/topaz/'+expid+'/'
+date = datetime(int(yyyy), int(mth), 1)
+
+lon_name='lon'
+lat_name='lat'
+if varname=='to':
+    zidx=1
+    grid_name=outdir+varname+'/'+fq+'/'+expid+'_'+varname+'_'+str(zidx)+'_'+fq+'_'+date.strftime('%Y%m')+'_MR25.nc'
+else:
+    grid_name=outdir+varname+'/'+fq+'/'+expid+'_'+varname+'_'+fq+'_'+date.strftime('%Y%m')+'_MR25.nc'
+
+wavelength = 700  #choice of spatial cutoff for high pass filter in km
+
+g = RegularGridDataset(grid_name, lon_name, lat_name) #assumes no time index, might only take the first time index
+g.bessel_high_filter(varname, wavelength, order=1) #perfroms only on 1 time index
+if varname=='to':
+    zidx=1
+    g.write(outdir+varname+'/'+fq+'/'+expid+'_'+varname+'_'+str(zidx)+'_'+fq+'_'+date.strftime('%Y%m')+'_MR25_hp'+str(wavelength)+'.nc')
+else:
+    g.write(outdir+varname+'/'+fq+'/'+expid+'_'+varname+'_'+fq+'_'+date.strftime('%Y%m')+'_MR25_hp'+str(wavelength)+'.nc')
+
+

eddy_composite/remap_mm_r2b8O_MR25.job

@@ -0,0 +1,71 @@
+#!/bin/bash
+#SBATCH --account=mh0256
+#SBATCH --job-name=remap25
+#SBATCH --partition=compute
+#SBATCH --nodes=1
+#SBATCH --mem=0
+#SBATCH --output=remap25.run.%j.o
+#SBATCH --error=remap25.run.%j.o
+#SBATCH --exclusive
+#SBATCH --time=08:00:00
+#SBATCH --mail-type=END
+#SBATCH --mail-user=dian.putrasahan@mpimet.mpg.de
+
+if [ -z "$1" ] ; then
+    echo "invalid number of parameters: need varname"
+    exit 1      #exit out of program
+fi
+if [ -z "$2" ] ; then
+    echo "invalid number of parameters: need expid"
+    exit 1      #exit out of program
+fi
+
+# ©2023 MPI-M, Dian Putrasahan
+
+varname=$1
+expid=$2
+
+module load cdo
+
+find_files=/home/m/m300466/pyfuncs/find_files
+
+#varname=zos
+#varname=Wind_Speed_10m
+#varname=atmos_fluxes_HeatFlux_Latent
+#varname=mlotst
+#varname=sea_level_pressure
+#varname=to
+#varname=atmos_fluxes_FrshFlux_Precipitation
+
+#expid=ngc2013
+#expid=rthk001
+if [ $expid == 'ngc2013' ]; then
+   catalog=/home/m/m300466/NGC2/intake/${expid}_bu1213.json
+elif [ $expid == 'rthk001' ]; then
+   catalog=/home/m/m300466/NGC2/intake/rthk001_test-sorted_stylized.json
+fi
+
+freqm=mm
+ocngrid=/pool/data/ICON/grids/public/mpim/0034/icon_grid_0034_R02B08_O.nc
+remapwgt=/work/mh0256/m300466/DPP/weights/r2b8O_MR25_yconremapweights.nc
+targetgrid=/work/mh0287/m300466/topaz/ngc2013/zos/zos_2054-01-01.nc
+
+
+#Obtain a list of path and file names that contain specified variable, frequency and time range (Feb 2020 - Dec 2049)
+varfilelist=$( $find_files --catalog_file=${catalog} ${varname} ${expid} --frequency=1month --time_range 2020-02-01 2050-02-01T23)
+echo $varfilelist | tr ' ' '\n' > filelist_0120
+
+#Extract upper 500m, monthly mean and remap
+mkdir -p /work/mh0287/m300466/topaz/${expid}/${varname}
+mkdir -p /work/mh0287/m300466/topaz/${expid}/${varname}/${freqm}
+
+if [ $varname == to ]; then
+   zidx=1
+   varfile=/work/mh0287/m300466/topaz/${expid}/${varname}/${freqm}/${expid}_${varname}_${zidx}_${freqm}_202002-205001_MR25.nc
+   cdo -L -P 128 -remap,${targetgrid},${remapwgt} -shifttime,-12hour -setgrid,${ocngrid} -select,name=${varname},levidx=${zidx} [ ${varfilelist} ] ${varfile}
+else
+   varfile=/work/mh0287/m300466/topaz/${expid}/${varname}/${freqm}/${expid}_${varname}_${freqm}_202002-205001_MR25.nc
+   cdo -L -P 128 -remap,${targetgrid},${remapwgt} -shifttime,-12hour -setgrid,${ocngrid} -select,name=${varname} [ ${varfilelist} ] ${varfile}
+fi
+#We shift the time stamp by -12hours in order to get it to the correct month
+

eddy_composite/split_mm_yrmth.sh

@@ -0,0 +1,93 @@
+#!/bin/bash
+# set -ex
+# ©2023 MPI-M, Dian Putrasahan
+
+#------------------------
+# users choice
+#------------------------
+
+#varname=zos
+#varname=Wind_Speed_10m
+#varname=atmos_fluxes_HeatFlux_Latent
+varname=mlotst
+#varname=to
+#varname=atmos_fluxes_FrshFlux_Precipitation
+
+expid=ngc2013
+#expid=rthk001
+
+freqm=mm
+zidx=1
+outdir=/work/mh0287/m300466/topaz/${expid}/${varname}/${freqm}
+
+#========================
+YYYYs=2020
+MMs=2
+YYYYe=2050
+MMe=1
+YYYYint=$YYYYs
+declare -a datearr
+declare -a dayarr
+declare -a yyarr
+declare -a mmarr
+#declare -a ddarr
+#dayarr=()
+
+while [ $YYYYint -le $YYYYe ]; do
+  if [ $YYYYint -eq $YYYYs ]; then
+     MMint=$MMs
+  else
+     MMint=1
+  fi
+  if [ $YYYYint -eq $YYYYe ]; then
+     MMend=$MMe
+  else
+     MMend=12
+  fi
+  while [ $MMint -le $MMend ]; do
+     if [ $MMint -le 9 ]; then
+        MMint=`eval echo "0$MMint"`
+     fi
+     #echo "Extract for yr$YYYYint mth$MMint day$DDint"
+     #dayarr+=( $fdate )
+     datearr+=( "${YYYYint}-${MMint}" )
+     dayarr+=( "${YYYYint}${MMint}" )
+     yyarr+=( "${YYYYint}" )
+     mmarr+=( "${MMint}" )
+     #ddarr+=( "${DDint}" )
+     MMint=`expr $MMint + 01`
+  done
+  YYYYint=`expr $YYYYint + 01`
+done
+
+#array length
+echo ${#dayarr[@]}
+##array list
+#echo ${dayarr[@]}
+
+for tidx in $(seq 0 ${#dayarr[@]});
+do
+  echo "ddate="${datearr[$tidx]}
+  ddate=${datearr[$tidx]}
+  echo "fdate="${dayarr[$tidx]}
+  fdate=${dayarr[$tidx]}
+  yyyy=${yyarr[$tidx]}
+  mm=${mmarr[$tidx]}
+  #Remove depth dimension in SST
+  varfile1=${outdir}/${varname}_${freqm}_${fdate}.nc
+  if [ ${varname} == 'to' ]; then
+     bigfile=${outdir}/${expid}_${varname}_${zidx}_${freqm}_202002-205001_MR25.nc
+     file2sm=${outdir}/${expid}_${varname}_${zidx}_${freqm}_${fdate}_MR25.nc
+     cdo -L -select,year=${yyyy},month=${mm} ${bigfile} ${varfile1}
+     ncwa -a depth ${varfile1} ${file2sm}
+     ncks -C -O -x -v depth ${file2sm} ${file2sm}
+     rm $varfile1
+  else
+     bigfile=${outdir}/${expid}_${varname}_${freqm}_202002-205001_MR25.nc
+     file2sm=${outdir}/${expid}_${varname}_${freqm}_${fdate}_MR25.nc
+     cdo -L -select,year=${yyyy},month=${mm} ${bigfile} ${file2sm}
+     #ncks -C -O -x -v lon,lat ${file2sm} ${file2sm}
+  fi
+done
+
+

eddy_composite/submit_IDeddy_mm.job

@@ -0,0 +1,39 @@
+#!/bin/sh
+#SBATCH --account=mh0256
+#SBATCH --job-name=IDeddy
+#SBATCH --partition=compute
+#SBATCH --nodes=1
+#SBATCH --mem=0
+#SBATCH --output=IDeddy_%j.log
+#SBATCH --error=IDeddy_%j.log
+#SBATCH --exclusive
+#SBATCH --time=02:00:00
+#SBATCH --mail-type=END
+#SBATCH --mail-user=dian.putrasahan@mpimet.mpg.de
+
+if [ -z "$1" ] ; then
+    echo "invalid number of parameters: need year"
+    exit 1      #exit out of program
+fi
+yr=$1
+
+module unload python3
+source /work/mh0256/m300466/miniconda3/bin/activate
+conda activate t2
+
+#for yr in $(seq 2050 2050);
+#for yr in $(seq 2020 2020);
+#for yr in $(seq 2021 2023);
+#for yr in $(seq 2024 2024);
+#for yr in $(seq 2025 2029);
+#for yr in $(seq 2030 2034);
+#for yr in $(seq 2035 2039);
+#for yr in $(seq 2040 2044);
+#for yr in $(seq 2045 2049);
+#do
+   for mth in $(seq 1 12);
+   #for mth in $(seq 1 1);
+   do
+      python IDeddy_mm.py ${yr} ${mth}
+   done
+#done