From ca3925fd5863f3be97d7f6aa793e67c94273492b Mon Sep 17 00:00:00 2001
From: voigta80 <aiko.voigt@univie.ac.at>
Date: Thu, 16 Sep 2021 14:42:09 +0200
Subject: [PATCH] Added benchmarking routines
---
benchmarks/benchmark_2d.py | 74 +++++++++++++++
benchmarks/benchmark_2d.run | 34 +++++++
benchmarks/benchmark_3d.py | 74 +++++++++++++++
benchmarks/benchmark_3d.run | 34 +++++++
benchmarks/benchmark_cubulation.ipynb | 126 --------------------------
benchmarks/benchmark_cubulation.py | 41 +++++++++
benchmarks/benchmark_cubulation.run | 36 ++++++++
benchmarks/benchmark_cubulation.txt | 14 +++
examples/find_radius.py | 25 +++++
examples/find_startcell.py | 35 +++++++
setup.py | 2 +-
11 files changed, 368 insertions(+), 127 deletions(-)
create mode 100644 benchmarks/benchmark_2d.py
create mode 100755 benchmarks/benchmark_2d.run
create mode 100644 benchmarks/benchmark_3d.py
create mode 100755 benchmarks/benchmark_3d.run
delete mode 100644 benchmarks/benchmark_cubulation.ipynb
create mode 100644 benchmarks/benchmark_cubulation.py
create mode 100755 benchmarks/benchmark_cubulation.run
create mode 100644 benchmarks/benchmark_cubulation.txt
create mode 100644 examples/find_radius.py
create mode 100644 examples/find_startcell.py
diff --git a/benchmarks/benchmark_2d.py b/benchmarks/benchmark_2d.py
new file mode 100644
index 0000000..7e45a2a
--- /dev/null
+++ b/benchmarks/benchmark_2d.py
@@ -0,0 +1,74 @@
+# Perform benchmarking of 2d cloud data given a previously computed cubulation
+
+# parse command line parameters
+import sys
+resol = sys.argv[1]
+startcell = int(sys.argv[2])
+searchrad = int(sys.argv[3])
+
+# gridfile including path
+gridfile = '/work/bb1018/b380459/NAWDEX/grids/icon-grid_nawdex_78w40e23n80n_'+resol+'.nc'
+
+# load other needed packages
+sys.path.append('/pf/b/b380459/connected-components-3d/')
+sys.path.append('/pf/b/b380459/tricco/')
+import tricco
+import datetime
+
+print(' ')
+print('-----------------------------------------')
+print('Working on resolution of', resol)
+
+# load previously computed cubulation
+import numpy as np
+cubulpath = '/scratch/b/b380459/tricco_output/'
+cubulfile = cubulpath+'/icon-grid_nawdex_78w40e23n80n_'+resol+'_cubulation_start'+str(startcell)+'_radius'+str(searchrad)+'.npy'
+cubulation = np.load(cubulfile, allow_pickle=True)
+
+# read in cloud data
+
+# cloud file depends on resolution
+datapath={'R80000m': '/work/bb1018/b380459/NAWDEX/ICON_OUTPUT_NWP/nawdexnwp-80km-mis-0001/',
+ 'R40000m': '/work/bb1018/b380459/NAWDEX/ICON_OUTPUT_NWP/nawdexnwp-40km-mis-0001/',
+ 'R20000m': '/work/bb1018/b380459/NAWDEX/ICON_OUTPUT_NWP/nawdexnwp-20km-mis-0001/',
+ 'R10000m': '/work/bb1018/b380459/NAWDEX/ICON_OUTPUT_NWP/nawdexnwp-10km-mis-0001/'}
+datafile={'R80000m': 'nawdexnwp-80km-mis-0001_2016092200_2d_30min_DOM01_ML_',
+ 'R40000m': 'nawdexnwp-40km-mis-0001_2016092200_2d_30min_DOM01_ML_',
+ 'R20000m': 'nawdexnwp-20km-mis-0001_2016092200_2d_30min_DOM01_ML_',
+ 'R10000m': 'nawdexnwp-10km-mis-0001_2016092200_2d_30min_DOM01_ML_'}
+
+dtime_dat = list()
+dtime_ver = list()
+dtime_edg = list()
+# loop over 1 day of 30-min output data --> 48 timesteps
+for time in range(10,59):
+ # read in data
+ begin_time = datetime.datetime.now()
+ field, field_cube = tricco.prepare_field(model='ICON', path=datapath[resol],
+ file=datafile[resol]+'00'+str(time)+'.nc',
+ var='clct', threshold=85.0, cubulation=cubulation)
+ end_time = datetime.datetime.now()
+ dtime_dat.append(end_time-begin_time)
+ # perform connected component analysis for vertex connectivity
+ begin_time = datetime.datetime.now()
+ _ = tricco.compute_connected_components_2d(cubulation, field_cube, connectivity = 'vertex')
+ end_time = datetime.datetime.now()
+ dtime_ver.append(end_time-begin_time)
+ # perform connected component analysis for edge connectivity
+ begin_time = datetime.datetime.now()
+ _ = tricco.compute_connected_components_2d(cubulation, field_cube, connectivity = 'edge')
+ end_time = datetime.datetime.now()
+ dtime_edg.append(end_time-begin_time)
+
+print(' ')
+print('Time to read in 1 time step of data, including moving it onto the 3d cubulated grid; done for 48 time steps')
+print('Mean:', np.mean(dtime_dat), 'Min:', np.min(dtime_dat), 'Max:', np.max(dtime_dat))
+print(' ')
+print('Time to do 2d connected component labeling with vertex connectivity for 1 time step; done for 48 time steps')
+print('Mean:', np.mean(dtime_ver), 'Min:', np.min(dtime_ver), 'Max:', np.max(dtime_ver))
+print(' ')
+print('Time to do 2d connected component labeling with edge connectivity for 1 time step: done for 48 time steps')
+print('Mean:', np.mean(dtime_edg), 'Min:', np.min(dtime_edg), 'Max:', np.max(dtime_edg))
+
+print('-----------------------------------------')
+print(' ')
diff --git a/benchmarks/benchmark_2d.run b/benchmarks/benchmark_2d.run
new file mode 100755
index 0000000..14d0cea
--- /dev/null
+++ b/benchmarks/benchmark_2d.run
@@ -0,0 +1,34 @@
+#!/bin/bash
+#=============================================================================
+# mistral batch job parameters
+#-----------------------------------------------------------------------------
+#SBATCH --account=bb1152
+#SBATCH --job-name=benchmark_2d.run
+#SBATCH --partition=compute
+#SBATCH --nodes=1
+#SBATCH --threads-per-core=1
+#SBATCH --output=/pf/b/b380459/BigDataClouds/tricco/benchmarks/LOG.benchmark_2d.run.%j.o
+#SBATCH --error=/pf/b/b380459/BigDataClouds/tricco/benchmarks/LOG.benchmark_2d.run.%j.o
+#SBATCH --exclusive
+#SBATCH --time=08:00:00
+
+cd /pf/b/b380459/BigDataClouds/tricco/benchmarks
+
+MYPYTHON="/pf/b/b380459/conda-envs/Nawdex-Hackathon/bin/python3.8"
+echo "Working with" ${MYPYTHON}
+
+for i in {1..10}; do
+ $MYPYTHON benchmark_2d.py R80000m 5738 102
+done
+
+for i in {1..10}; do
+ $MYPYTHON benchmark_2d.py R40000m 18538 230
+done
+
+for i in {1..10}; do
+ $MYPYTHON benchmark_2d.py R20000m 69309 460
+done
+
+for i in {1..10}; do
+ $MYPYTHON benchmark_2d.py R10000m 264792 2000
+done
diff --git a/benchmarks/benchmark_3d.py b/benchmarks/benchmark_3d.py
new file mode 100644
index 0000000..7687eb8
--- /dev/null
+++ b/benchmarks/benchmark_3d.py
@@ -0,0 +1,74 @@
+# Perform benchmarking of 3d cloud data given a previously computed cubulation
+
+# parse command line parameters
+import sys
+resol = sys.argv[1]
+startcell = int(sys.argv[2])
+searchrad = int(sys.argv[3])
+
+# gridfile including path
+gridfile = '/work/bb1018/b380459/NAWDEX/grids/icon-grid_nawdex_78w40e23n80n_'+resol+'.nc'
+
+# load other needed packages
+sys.path.append('/pf/b/b380459/connected-components-3d/')
+sys.path.append('/pf/b/b380459/tricco/')
+import tricco
+import datetime
+
+print(' ')
+print('-----------------------------------------')
+print('Working on resolution of', resol)
+
+# load previously computed cubulation
+import numpy as np
+cubulpath = '/scratch/b/b380459/tricco_output/'
+cubulfile = cubulpath+'/icon-grid_nawdex_78w40e23n80n_'+resol+'_cubulation_start'+str(startcell)+'_radius'+str(searchrad)+'.npy'
+cubulation = np.load(cubulfile, allow_pickle=True)
+
+# read in cloud data
+
+# cloud file depends on resolution
+datapath={'R80000m': '/work/bb1018/b380459/NAWDEX/ICON_OUTPUT_NWP/nawdexnwp-80km-mis-0001/',
+ 'R40000m': '/work/bb1018/b380459/NAWDEX/ICON_OUTPUT_NWP/nawdexnwp-40km-mis-0001/',
+ 'R20000m': '/work/bb1018/b380459/NAWDEX/ICON_OUTPUT_NWP/nawdexnwp-20km-mis-0001/',
+ 'R10000m': '/work/bb1018/b380459/NAWDEX/ICON_OUTPUT_NWP/nawdexnwp-10km-mis-0001/'}
+datafile={'R80000m': 'nawdexnwp-80km-mis-0001_2016092200_3dcloud_DOM01_ML_',
+ 'R40000m': 'nawdexnwp-40km-mis-0001_2016092200_3dcloud_DOM01_ML_',
+ 'R20000m': 'nawdexnwp-20km-mis-0001_2016092200_3dcloud_DOM01_ML_',
+ 'R10000m': 'nawdexnwp-10km-mis-0001_2016092200_3dcloud_DOM01_ML_'}
+
+dtime_dat = list()
+dtime_ver = list()
+dtime_edg = list()
+# loop over 1 day of 30-min output data --> 48 timesteps
+for time in range(10,59):
+ # read in data
+ begin_time = datetime.datetime.now()
+ field, field_cube = tricco.prepare_field_lev(model='ICON', path=datapath[resol],
+ file=datafile[resol]+'00'+str(time)+'.nc',
+ var='clc', threshold=85.0, cubulation=cubulation)
+ end_time = datetime.datetime.now()
+ dtime_dat.append(end_time-begin_time)
+ # perform connected component analysis for vertex connectivity
+ begin_time = datetime.datetime.now()
+ _ = tricco.compute_connected_components_3d(cubulation, field_cube, connectivity = 'vertex')
+ end_time = datetime.datetime.now()
+ dtime_ver.append(end_time-begin_time)
+ # perform connected component analysis for edge connectivity
+ begin_time = datetime.datetime.now()
+ _ = tricco.compute_connected_components_3d(cubulation, field_cube, connectivity = 'edge')
+ end_time = datetime.datetime.now()
+ dtime_edg.append(end_time-begin_time)
+
+print(' ')
+print('Time to read in 1 time step of data, including moving it onto the 3d x lev cubulated grid; done for 48 time steps')
+print('Mean:', np.mean(dtime_dat), 'Min:', np.min(dtime_dat), 'Max:', np.max(dtime_dat))
+print(' ')
+print('Time to do 3d connected component labeling with vertex connectivity for 1 time step; done for 48 time steps')
+print('Mean:', np.mean(dtime_ver), 'Min:', np.min(dtime_ver), 'Max:', np.max(dtime_ver))
+print(' ')
+print('Time to do 3d connected component labeling with edge connectivity for 1 time step; done for 48 time steps')
+print('Mean:', np.mean(dtime_edg), 'Min:', np.min(dtime_edg), 'Max:', np.max(dtime_edg))
+
+print('-----------------------------------------')
+print(' ')
diff --git a/benchmarks/benchmark_3d.run b/benchmarks/benchmark_3d.run
new file mode 100755
index 0000000..85bfc33
--- /dev/null
+++ b/benchmarks/benchmark_3d.run
@@ -0,0 +1,34 @@
+#!/bin/bash
+#=============================================================================
+# mistral batch job parameters
+#-----------------------------------------------------------------------------
+#SBATCH --account=bb1152
+#SBATCH --job-name=benchmark_3d.run
+#SBATCH --partition=compute
+#SBATCH --nodes=1
+#SBATCH --threads-per-core=1
+#SBATCH --output=/pf/b/b380459/BigDataClouds/tricco/benchmarks/LOG.benchmark_3d.run.%j.o
+#SBATCH --error=/pf/b/b380459/BigDataClouds/tricco/benchmarks/LOG.benchmark_3d.run.%j.o
+#SBATCH --exclusive
+#SBATCH --time=08:00:00
+
+cd /pf/b/b380459/BigDataClouds/tricco/benchmarks
+
+MYPYTHON="/pf/b/b380459/conda-envs/Nawdex-Hackathon/bin/python3.8"
+echo "Working with" ${MYPYTHON}
+
+#for i in {1..10}; do
+# $MYPYTHON benchmark_3d.py R80000m 5738 102
+#done
+
+#for i in {1..10}; do
+# $MYPYTHON benchmark_3d.py R40000m 18538 230
+#done
+
+#for i in {1..10}; do
+# $MYPYTHON benchmark_3d.py R20000m 69309 460
+#done
+
+for i in {1..10}; do
+ $MYPYTHON benchmark_3d.py R10000m 264792 2000
+done
diff --git a/benchmarks/benchmark_cubulation.ipynb b/benchmarks/benchmark_cubulation.ipynb
deleted file mode 100644
index 0506e40..0000000
--- a/benchmarks/benchmark_cubulation.ipynb
+++ /dev/null
@@ -1,126 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Notebook for benchmarking runtime of cubulation\n",
- "\n",
- "Tests are performed on an exclusive compute node of the DKRZ supercomputer Mistral in Hamburg, Germany.\n",
- "\n",
- "A compute node has the following specs (https://www.dkrz.de/up/systems/mistral/configuration):\n",
- " * 2x 12-core Intel Xeon E5-2680 v3 (Haswell) @ 2.5GHz\n",
- " * 24 cores (48 logical CPUs)\n",
- " * 64 GB main memory\n",
- " \n",
- "As for the grids we consider limited-area ICON grids that cover a large part of the North Atlantic. They were for example used for the NAWDEX simulations described Senf, F., A. Voigt et al, 2020: Increasing Resolution and Resolving Convection Improve the Simulation of Cloudâ€Radiative Effects Over the North Atlantic, JGR Atmospheres. https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020JD032667."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Load required packages. Note: adding cc3d to system path is needed because it is required by tricco."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 1,
- "metadata": {},
- "outputs": [],
- "source": [
- "import timeit\n",
- "import sys\n",
- "sys.path.append('/pf/b/b380459/connected-components-3d/')\n",
- "sys.path.append('/pf/b/b380459/BigDataClouds/tricco/')\n",
- "import tricco"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Define the start triangle and radius of outward search as a function of grid resolution."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "metadata": {},
- "outputs": [],
- "source": [
- "dict_start = {'80000m': 5783, '40000m': 18538, '20000m': 69309, '10000m': 264792}\n",
- "dict_radius = {'80000m': 102 , '40000m': 230 , '20000m': 460 , '10000m': 2000 }"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Working on resolution of 80000m\n",
- "Time to read the grid:\n",
- "310 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n",
- "Time to compute cubulation:\n",
- "11.7 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n",
- "Working on resolution of 40000m\n",
- "Time to read the grid:\n",
- "163 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n",
- "Time to compute cubulation:\n",
- "1min 13s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n",
- "Working on resolution of 20000m\n",
- "Time to read the grid:\n",
- "144 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n",
- "Time to compute cubulation:\n",
- "11min 34s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n",
- "Working on resolution of 10000m\n",
- "Time to read the grid:\n",
- "312 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n",
- "Time to compute cubulation:\n"
- ]
- }
- ],
- "source": [
- "# path to grid files\n",
- "gridpath = '/work/bb1018/b380459/NAWDEX/grids/'\n",
- "\n",
- "# model resolutions under investigation\n",
- "for res in ['80000m', '40000m', '20000m', '10000m']:\n",
- " \n",
- " print('Working on resolution of', res)\n",
- " \n",
- " gridfile = 'icon-grid_nawdex_78w40e23n80n_R'+res+'.nc'\n",
- " print('Time to read the grid:')\n",
- " %timeit -r 1 -n 1 tricco.grid_functions.grid = tricco.prepare_grid(model='ICON',path=gridpath, file=gridfile)\n",
- " \n",
- " print('Time to compute cubulation:')\n",
- " %timeit -r 1 -n 1 tricco.compute_cubulation(start_triangle=dict_start[res], radius=dict_radius[res], print_progress=False)"
- ]
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Nawdex-Hackathon",
- "language": "python",
- "name": "nawdex-hackathon"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.8.10"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
diff --git a/benchmarks/benchmark_cubulation.py b/benchmarks/benchmark_cubulation.py
new file mode 100644
index 0000000..68bac40
--- /dev/null
+++ b/benchmarks/benchmark_cubulation.py
@@ -0,0 +1,41 @@
+# Perform benchmarking of cubulation routine on DKRZ Mistral
+
+# parse command line parameters
+import sys
+resol = sys.argv[1]
+startcell = int(sys.argv[2])
+searchrad = int(sys.argv[3])
+savecubul = int(sys.argv[4])
+
+# gridfile including path
+gridfile = '/work/bb1018/b380459/NAWDEX/grids/icon-grid_nawdex_78w40e23n80n_'+resol+'.nc'
+
+# load other needed packages
+sys.path.append('/pf/b/b380459/connected-components-3d/')
+sys.path.append('/pf/b/b380459/tricco/')
+import tricco
+import datetime
+
+print(' ')
+print('-----------------------------------------')
+print('Working on resolution of', resol)
+print('Start cell:', startcell, 'Search radius:', searchrad, 'Save cubulation:', savecubul)
+
+begin_time = datetime.datetime.now()
+tricco.grid_functions.grid = tricco.prepare_grid(model='ICON',path='/', file=gridfile)
+end_time = datetime.datetime.now()
+print('Time to read the grid:', end_time-begin_time)
+
+begin_time = datetime.datetime.now()
+cubulation = tricco.compute_cubulation(start_triangle=startcell, radius=searchrad, print_progress=False)
+end_time = datetime.datetime.now()
+print('Time to compute cubulation:', end_time-begin_time)
+
+# optional: save cubulation
+if savecubul == 1:
+ import numpy as np
+ np.save('/scratch/b/b380459/tricco_output/icon-grid_nawdex_78w40e23n80n_'
+ +resol+'_cubulation_start'+str(startcell)+'_radius'+str(searchrad), cubulation)
+
+print('-----------------------------------------')
+print(' ')
diff --git a/benchmarks/benchmark_cubulation.run b/benchmarks/benchmark_cubulation.run
new file mode 100755
index 0000000..9e0da59
--- /dev/null
+++ b/benchmarks/benchmark_cubulation.run
@@ -0,0 +1,36 @@
+#!/bin/bash
+#=============================================================================
+# mistral batch job parameters
+#-----------------------------------------------------------------------------
+#SBATCH --account=bb1152
+#SBATCH --job-name=benchmark_cubulation.run
+#SBATCH --partition=compute
+#SBATCH --nodes=1
+#SBATCH --threads-per-core=1
+#SBATCH --output=/pf/b/b380459/BigDataClouds/tricco/benchmarks/LOG.benchmark_cubulation.run.%j.o
+#SBATCH --error=/pf/b/b380459/BigDataClouds/tricco/benchmarks/LOG.benchmark_cubulation.run.%j.o
+#SBATCH --exclusive
+#SBATCH --time=08:00:00
+
+cd /pf/b/b380459/BigDataClouds/tricco/benchmarks
+
+MYPYTHON="/pf/b/b380459/conda-envs/Nawdex-Hackathon/bin/python3.8"
+echo "Working with" ${MYPYTHON}
+
+for i in {1..10}; do
+# $MYPYTHON benchmark_cubulation.py R80000m 5738 102 1
+ $MYPYTHON benchmark_cubulation.py R80000m 5568 200 1
+done
+
+for i in {1..10}; do
+# $MYPYTHON benchmark_cubulation.py R40000m 18538 230 1
+ $MYPYTHON benchmark_cubulation.py R40000m 18493 400 1
+done
+
+#for i in {1..10}; do
+# $MYPYTHON benchmark_cubulation.py R20000m 69309 460 1
+#done
+
+#for i in {1..10}; do
+# $MYPYTHON benchmark_cubulation.py R10000m 264792 2000 1
+#done
diff --git a/benchmarks/benchmark_cubulation.txt b/benchmarks/benchmark_cubulation.txt
new file mode 100644
index 0000000..737bdad
--- /dev/null
+++ b/benchmarks/benchmark_cubulation.txt
@@ -0,0 +1,14 @@
+Benchmarking runtime of cubulation
+
+Tests are performed on an exclusive compute node of the DKRZ supercomputer Mistral in Hamburg, Germany.
+
+A compute node has the following specs (https://www.dkrz.de/up/systems/mistral/configuration):
+ * 2x12-core Intel Xeon E5-2680 v3 (Haswell) @ 2.5GHz,
+ * 24 cores (48 logical CPUs),
+ * 64 GB main memory.
+
+As for the grids we consider limited-area ICON grids that cover a large part of the North Atlantic. They were for example used for the NAWDEX simulations described Senf, F., A. Voigt et al, 2020: Increasing Resolution and Resolving Convection Improve the Simulation of Cloudâ€Radiative Effects Over the North Atlantic, JGR Atmospheres. https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020JD032667.
+
+Usage:
+
+benchmark_cubulation.run is a batch job that is submitted to an exclusive compute node via sbatch and that calls benchmark_cubulation.py. The grid resolution, start triangle and search radius for the cubulation are handed over to benchmark_cubulation.py, which reads in the grid file and does the cubulation.
diff --git a/examples/find_radius.py b/examples/find_radius.py
new file mode 100644
index 0000000..08bd7c6
--- /dev/null
+++ b/examples/find_radius.py
@@ -0,0 +1,25 @@
+# Purpose: For a given start cell find the smallest radius that allows one to cover all grid cell.
+# This is helpful as a too large radius artifically inflates the size of the cubulated grid.
+
+# Written for the limited-area grids of ICON used in the Tricco introduction paper.
+
+# parse command line parameters
+import sys
+resol = sys.argv[1]
+start = int(sys.argv[2])
+radius = int(sys.argv[3])
+
+print('-----------------------------------------')
+print('Working on ICON grid with resolution', resol)
+print('Start cell :', start)
+print('Search radius:', radius)
+
+import sys
+sys.path.append('/pf/b/b380459/connected-components-3d/')
+sys.path.append('/pf/b/b380459/BigDataClouds/tricco/')
+import tricco
+
+tricco.grid_functions.grid = tricco.prepare_grid(model='ICON', path='./data/',
+ file='icon-grid_nawdex_78w40e23n80n_'+resol+'.nc')
+
+cubulation = tricco.compute_cubulation(start_triangle=start, radius=radius, print_progress=True)
diff --git a/examples/find_startcell.py b/examples/find_startcell.py
new file mode 100644
index 0000000..f27787d
--- /dev/null
+++ b/examples/find_startcell.py
@@ -0,0 +1,35 @@
+# Purpose: Find the cell closest to a given latitude-longitude position.
+# This is helpful to set the start cell of the cubulation routine.
+
+# Written for the limited-area grids of ICON used in the Tricco introduction paper.
+
+# convert rad to deg
+import numpy as np
+rad2deg=180.0/np.pi
+
+# parse command line parameters
+import sys
+resol = sys.argv[1]
+lat = float(sys.argv[2])
+lon = float(sys.argv[3])
+
+print('-----------------------------------------')
+print('Working on ICON grid with resolution', resol)
+print('Searching for cell closest to lat', lat, 'and lon', lon)
+
+# gridfile including path
+gridfile = '/work/bb1018/b380459/NAWDEX/grids/icon-grid_nawdex_78w40e23n80n_'+resol+'.nc'
+
+# load lat-lon info of grid and convert to deg
+import xarray as xr
+ds_grid = xr.load_dataset(gridfile)
+clat = rad2deg*ds_grid['clat'].values
+clon = rad2deg*ds_grid['clon'].values
+
+dist = np.power(clat-lat,2) + np.power(clon-lon,2)
+
+print('Closest cell has index', np.argmin(dist))
+print('Note: The startcell for tricco is the cell index - 1.')
+print(' This is because the found cell index is on the ICON grid')
+print(' and the ICON indexing starts with 1.')
+print('-----------------------------------------')
diff --git a/setup.py b/setup.py
index b17f04e..b69b503 100644
--- a/setup.py
+++ b/setup.py
@@ -1,6 +1,6 @@
from setuptools import setup, find_packages
-VERSION = '0.0.4'
+VERSION = '0.0.1'
DESCRIPTION = 'TriCCo'
LONG_DESCRIPTION = 'TriCCo: a python package for connected component labeling on triangular grids'
--
GitLab