Repository: jhardenberg/smmregrid (language: python)

Issue:
Refactoring of smmregrid to support for more complex dataset structures
There is the growing evidence that smmregrid is not flexible enough to support to ongoing stress coming from external application. This issue will serve as a landing point for the develpoment that we are planning to bring in order to make it more flexible and sustainable. 

Most important limitation that I am aware of is the impossibility of dealing with xarray dataset which does not share dimensions in exact way. You cannot remap a source which have both oceanic and atmospheric data. 

In order to overcome this limitation, the idea is make smmregrid work based on multiple `gridtype` instead of assuming that the data is all on the same grids. This is somehow not different from what CDO is doing.

1. Introduction of a function to identify the `gridtype` that are available in a specific Dataset/DataArray based on shared dimensions, and exclude those which are not relevant (the ones that include `bnds` for example. Please note that this might be required by CDO so we should learn how to bring them along). Any `gridtype` can be identifed based on a tuple of dimensions (no strict naming required). Then, we can build a dictionary and exploit all the tools we have already available to detect vertical and horizontal dimensions, as well as variables which are lying on the specific `gridtype`.  This is planned to be a class and I did already some successful test offline.
2. Make all the smmregrid class object based on  the `gridtype`. So the `cdo_generate_weights` will return not a simple xarray object, but a dictionary associated with each `gridtype`. All the object of the class are now expected to work with dictionaries, so that we can point to the required grid everytime.
3. As a final step, the `regrid` call should check what is the `gridtype` of the data that is fed into the method, and match it with the ones available in weights. 

Overall, this should massively improve the flexibility of the tool as well as the handling of the vertical coordinates, which will come now as inner property of the `gridtype`.
I plan also to move the `apply_weights` into the `Regrid` class to minimize the redundancy of the variables.

Tricky issues that I see:
- Deal with precomputed weights stored on disk in the correct way, so that we associate each of them to the right `gridtype`. This is probably something more from AQUA but we should be think about it. 
- Do not mess with DataArray and Dataset: currently the core regrid functions works on DataArray but having multiple grids must require the proper handling of the Dataset
- Find the right position where to load the data: we need to access metadata quite early so even if file are supplied this has to be moved at the init of the class.
- More that I do not see now

I did a test in my local branch and it seems feasible to achieve. I already have a nice `GridInspector` class that provides most of the grid information. Still, I very far from having something decent so I am not opening a PR yet. 
Speed evaluation of the smmregrid tool
This issue is to keep track of the speed tests that I have been to see what is the optimal configuration for the regridder based on #1 

The tests are based on files on different grids (curvilinear, gaussian, gaussian reduced, lonlat and unstructured) to cover all the possibilities, with 2D files, files with mask (i.e. ocean files) and files with pressure levels. We also tested the access of the entire xarray.Dataset versus working on the single xarray.DataArray. The writing of the NetCDF output is also assessed. All tests are run with conservative remapping.

The tests can be found in the playground notebook, and are based on multiple repetition (usually 20 fo each operation). https://github.com/jhardenberg/smmregrid/blob/devel/extend/playground.ipynb 


Gridtype revisit
I found a tricky bug which I am trying to work on here. The problem is that GridType object have an equality based on dimensions, horizontal+vertical+extra. However, we do not need this while we should always consider that horizontal and vertical case only. 

However, the issue is that variable identification should not a feature of GridType but rather of GridInspector. I am trying to revisit this here but it will require a bit of work. I assume we could merge the current refactoring and then further dig with this. 


Discussion / hints:

Commit b08b045f6fbb6cf63d439676b207e702d182181a establish a good starting point:

|                 |   NVars | NRecords          |   CDO |   SMM (Dataset) |   SMM (DataArray) |   SMM (DataArray+NoMask) |   SMM (Dataset+Write) |   SMM (DataArray+Write) |
|:----------------|--------:|:------------------|------:|----------------:|------------------:|-------------------------:|----------------------:|------------------------:|
| onlytos-ipsl.nc |       1 | (12, 332, 362)    |     1 |        0.216799 |       0.0778903   |              0.0504366   |              0.997036 |                0.77708  |
| tas-ecearth.nc  |       1 | (12, 256, 512)    |     1 |        0.226347 |       0.0900398   |              0.0631857   |              1.14144  |                0.958676 |
| 2t-era5.nc      |       1 | (12, 73, 144)     |     1 |        0.170659 |       0.094557    |              0.0610341   |              0.845712 |                0.765937 |
| tos-fesom.nc    |       1 | (12, 126859)      |     1 |        0.113976 |       0.0399258   |              0.0256824   |              0.755877 |                0.623671 |
| ua-ecearth.nc   |       1 | (2, 19, 256, 512) |     1 |        0.398825 |       0.0677817   |              0.0441846   |              1.61761  |                1.35922  |
| mix-cesm.nc     |       4 | (12, 192, 288)    |     1 |        0.549228 |       0.0688021   |              0.0452126   |              1.73817  |                0.670783 |
| era5-mon.nc     |       1 | (864, 721, 1440)  |     1 |        0.825034 |       0.000651573 |              0.000439713 |              1.883    |                1.0852   |

Few points:

- There is still too much overhead in the Dataset call. This is likely due to some safecheck we are using for the masks. 
- Masks count for about 50% more computational increase. So far it is unsafe to run without masks.
- Performance are very good for big files (see the last era5-montly interpolation) which suggest great scaling
- Writing files bring we back to square one, putting on the same line as CDO (but this is expected).

**IMPORTANT: this numbers does not take into account the loading of the data**
Numbers are less incredible if we take into account the loading of the data into memory i.e. `xarray.load()`

|                 |   NVars | NRecords          |   CDO |   CDO (NoLoad) |   SMM (Dataset) |   SMM (DataArray) |   SMM (DataArray+NoLoad) |   SMM (DataArray+NoMask) |   SMM (Dataset+Write) |   SMM (DataArray+Write) |
|:----------------|--------:|:------------------|------:|---------------:|----------------:|------------------:|-------------------------:|-------------------------:|----------------------:|------------------------:|
| onlytos-ipsl.nc |       1 | (12, 332, 362)    |     1 |       0.85244  |        0.468714 |          0.36703  |                0.0898978 |                 0.261572 |              0.451468 |                0.436607 |
| tas-ecearth.nc  |       1 | (12, 256, 512)    |     1 |       0.976596 |        0.556708 |          0.498547 |                0.103215  |                 0.35501  |              0.583348 |                0.561232 |
| 2t-era5.nc      |       1 | (12, 73, 144)     |     1 |       1.00737  |        0.295552 |          0.250174 |                0.0685831 |                 0.163858 |              0.354942 |                0.341436 |
| tos-fesom.nc    |       1 | (12, 126859)      |     1 |       0.975193 |        0.354071 |          0.321254 |                0.052417  |                 0.226885 |              0.370075 |                0.367297 |
| ua-ecearth.nc   |       1 | (2, 19, 256, 512) |     1 |       0.886853 |        0.623172 |          0.597676 |                0.171645  |                 0.476483 |              0.834928 |                0.849127 |
| mix-cesm.nc     |       4 | (12, 192, 288)    |     1 |       0.828504 |        0.903087 |          0.246191 |                0.0574815 |                 0.172714 |              1.21837  |                0.3044   |
| era5-mon.nc     |       1 | (864, 721, 1440)  |     1 |       1.00124  |        0.68279  |          0.688481 |                0.338506  |                 0.66355  |              0.696658 |                0.705815 |

We are still faster than CDO for single DataArray, but the speedup is small. 
Conversely, we get to very bad scaling when we use `dask`. This is pretty unclear why.

|                 |   Workers |     CDO |   Dask Compute |   Dask Load |
|:----------------|----------:|--------:|---------------:|------------:|
| onlytos-ipsl.nc |         0 | 1.33525 |       0.802709 |    0.80034  |
| tas-ecearth.nc  |         0 | 1.31848 |       1.04944  |    0.976837 |
| onlytos-ipsl.nc |         1 | 1.36808 |       3.72327  |    2.29133  |
| tas-ecearth.nc  |         1 | 1.36115 |       2.75875  |    2.80817  |
| onlytos-ipsl.nc |         2 | 1.3585  |       4.1103   |    2.21112  |
| tas-ecearth.nc  |         2 | 1.40967 |       2.85582  |    2.78464  |
| onlytos-ipsl.nc |         8 | 1.356   |       7.77573  |    5.73662  |
| tas-ecearth.nc  |         8 | 1.43555 |       3.41429  |    2.84304  |
Last commit in #2 suggest significant improvements. Considering that we are not using dask yet, this can be considered as a success. 

|                 |   NVars | NRecords          |   CDO |   CDO (NoLoad) |   SMM (Dataset) |   SMM (DataArray) |   SMM (DataArray+NoLoad) |   SMM (Dataset+Write) |   SMM (DataArray+Write) |
|:----------------|--------:|:------------------|------:|---------------:|----------------:|------------------:|-------------------------:|----------------------:|------------------------:|
| onlytos-ipsl.nc |       1 | (12, 332, 362)    |     1 |       0.924445 |        0.377513 |          0.347206 |                0.078164  |              0.429603 |                0.431012 |
| tas-ecearth.nc  |       1 | (12, 256, 512)    |     1 |       0.927699 |        0.410739 |          0.38098  |                0.0854462 |              0.468425 |                0.461594 |
| 2t-era5.nc      |       1 | (12, 73, 144)     |     1 |       0.924335 |        0.204548 |          0.160671 |                0.0439383 |              0.253002 |                0.248918 |
| tos-fesom.nc    |       1 | (12, 126859)      |     1 |       1.01328  |        0.332096 |          0.32686  |                0.048309  |              0.379621 |                0.372826 |
| ua-ecearth.nc   |       1 | (2, 19, 256, 512) |     1 |       0.898331 |        0.508237 |          0.475537 |                0.156023  |              0.74989  |                0.735353 |
| mix-cesm.nc     |       4 | (12, 192, 288)    |     1 |       0.886593 |        0.652437 |          0.181319 |                0.0456369 |              0.92309  |                0.248312 |
| era5-mon.nc     |       1 | (864, 721, 1440)  |     1 |       0.990528 |        0.757665 |          0.794592 |                0.33693   |              0.822056 |                0.803797 |


Fix patch:
diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
index bdaab28..b7a704b 100644
--- a/.github/workflows/python-publish.yml
+++ b/.github/workflows/python-publish.yml
@@ -21,7 +21,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - name: Set up Python
       uses: actions/setup-python@v3
       with:
diff --git a/.gitignore b/.gitignore
index 3dc8f23..f0dc778 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,5 @@ smmregrid.egg-info
 dist
 __pycache__
 *.idx
+docs/sphinx/build
+.coverage*
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
new file mode 100644
index 0000000..848a79a
--- /dev/null
+++ b/.readthedocs.yaml
@@ -0,0 +1,28 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "mambaforge-4.10"
+
+# set up the environment
+conda:
+  environment: environment.yml
+
+# install the local package
+python:
+  install:
+    - method: pip
+      path: .
+      extra_requirements:
+        - docs
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+   configuration: docs/sphinx/source/conf.py
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2596eec..a2d7a00 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,16 @@
 All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 
+## [v0.1.0]
+
+- Refactoring of the grid handling to possibly support more complex data structures (#33)
+- Introduction of `GridType()` class to handle possible multiple grids in the future (#33)
+- Introduction of `GridInspector()` class to investigate properties of a Dataset/DataArray (#33)
+- Minimal documentation via Sphinx available through ReadTheDocs (#33)
+- Preliminar support to target unstructured grids (#33)
+- `CdoGenerate()` class to support weights generator through `weights()` method (#33)
+- `CdoGenerate().areas()` method to generate areas based on CDO (#33)
+
 ## [v0.0.7]
 
 - Allow for `cdo_options` and move from `extra` to `cdo_extra` (#31)
diff --git a/README.md b/README.md
index 5866ffa..e6c592b 100644
--- a/README.md
+++ b/README.md
@@ -2,6 +2,8 @@
 [![PyTest](https://github.com/jhardenberg/smmregrid/actions/workflows/mambatest.yml/badge.svg)](https://github.com/jhardenberg/smmregrid/actions/workflows/mambatest.yml)
 [![Coverage Status](https://coveralls.io/repos/github/jhardenberg/smmregrid/badge.svg?branch=main)](https://coveralls.io/github/jhardenberg/smmregrid?branch=main)
 [![PyPI version](https://badge.fury.io/py/smmregrid.svg)](https://badge.fury.io/py/smmregrid)
+[![Documentation](https://readthedocs.org/projects/smmregrid/badge/?version=latest)](https://smmregrid.readthedocs.io/en/latest/?badge=latest)
+      
 
 # smmregrid
 A compact regridder using sparse matrix multiplication
diff --git a/docs/sphinx/Makefile b/docs/sphinx/Makefile
new file mode 100644
index 0000000..69fe55e
--- /dev/null
+++ b/docs/sphinx/Makefile
@@ -0,0 +1,19 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = sphinx-build
+SOURCEDIR     = source
+BUILDDIR      = build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
new file mode 100644
index 0000000..92ca1ce
--- /dev/null
+++ b/docs/sphinx/source/api.rst
@@ -0,0 +1,9 @@
+API Reference
+=============
+
+This section provides a detailed reference for smmregrid's Application Programming Interface (API).
+
+.. automodule:: smmregrid
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/docs/sphinx/source/conf.py b/docs/sphinx/source/conf.py
new file mode 100644
index 0000000..7147225
--- /dev/null
+++ b/docs/sphinx/source/conf.py
@@ -0,0 +1,45 @@
+# -*- coding: utf-8 -*-
+#
+# Configuration file for the Sphinx documentation builder.
+#
+# This file does only contain a selection of the most common options. For a
+# full list see the documentation:
+# http://www.sphinx-doc.org/en/master/config
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+import os
+import sys
+sys.path.insert(0, os.path.abspath('../../../'))
+
+# -- Project information -----------------------------------------------------
+
+project = 'smmregrid'
+copyright = '2022, Jost von Hardenberg and Paolo Davini'
+author = 'Jost von Hardenberg and Paolo Davini'
+
+# -- General configuration ---------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
+
+extensions = ['sphinx.ext.autodoc', 'sphinx.ext.napoleon']
+napoleon_google_docstring = True
+napoleon_numpy_docstring = False
+
+templates_path = ["_templates"]
+autoclass_content='both'
+exclude_patterns = []
+
+
+# -- Options for HTML output -------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
+
+html_theme = "sphinx_rtd_theme"
+# html_static_path = ["_static"]
+html_theme_options = {
+    "collapse_navigation": False,
+    "sticky_navigation": True,
+    "navigation_depth": 4,
+}
\ No newline at end of file
diff --git a/docs/sphinx/source/index.rst b/docs/sphinx/source/index.rst
new file mode 100644
index 0000000..fb25fec
--- /dev/null
+++ b/docs/sphinx/source/index.rst
@@ -0,0 +1,15 @@
+.. ECmean4 documentation master file, created by
+   sphinx-quickstart on Tue May  3 10:58:01 2022.
+   You can adapt this file completely to your liking, but it should at least
+   contain the root `toctree` directive.
+
+Welcome to smmregrid's documentation!
+=====================================
+
+.. toctree::
+   :maxdepth: 2
+
+   intro
+   installation
+   usage
+   api
\ No newline at end of file
diff --git a/docs/sphinx/source/installation.rst b/docs/sphinx/source/installation.rst
new file mode 100644
index 0000000..0b8b29e
--- /dev/null
+++ b/docs/sphinx/source/installation.rst
@@ -0,0 +1,67 @@
+Installation
+============
+
+
+smmregrid is a lightweight python package, but it depends on ``cdo`` for weights installation, thus both installation options requires conda/mamba. 
+We recommend to use `mamba <https://mamba.readthedocs.io/en/latest/user_guide/mamba.html>`_ since it provides a lighter and deal in a better way with dependencies.
+
+Using PyPi
+----------
+
+.. warning::
+
+    Please note that although smmregrid is distributed via PyPi, it depends on packages that currently are available 
+    only on conda-forge and on configuration files available from the GitHub repository.
+    Therefore, the installation via pip requires the creation of a conda environment as well as the clone from the repository.
+
+
+It will bring you the last version available on PyPi.
+You can create a conda/mamba environment which incudes the python, `eccodes <https://github.com/ecmwf/eccodes-python>`_ and `cdo <https://code.mpimet.mpg.de/projects/cdo/>`_ dependencies, and then install smmregrid.
+However, you should start by cloning the repository from GitHub, since the configuration files used for running ECmean4 are placed there ::
+
+    > git clone https://github.com/jvonhard/smmregrid.git
+    > mamba create -n smmregrid "python>=3.8" cdo eccodes
+    > mamba activate smmregrid
+    > pip install smmregrid
+
+
+Using GitHub
+------------
+
+This method will allow you to have access at the most recent ECmean4 version but it requires a bit more of effort.
+
+As before, should clone from the Github Repository ::
+
+    > git clone https://github.com/jvonhard/smmregrid.git
+    
+.. note ::
+
+    Please note that if you clone with HTTPS you will not be able to contribute to the code, even if you are listed as collaborator.
+    If you want to be a developer you should clone with SSH and you should add your own SSH key on the GitHub portal: 
+    please check the `procedure on the Github website <https://docs.github.com/en/authentication/connecting-to-github-with-ssh/adding-a-new-ssh-key-to-your-github-account>`_ .
+
+Then you can through the ECmean4 folder ::
+
+    > cd smmregrid
+
+and then you can set up the conda/mamba environment ::
+
+    > mamba env create --name smmregrid -f environment.yml
+
+Then you should activate the environment ::
+
+    > mamba activate smmregrid
+
+
+Requirements
+------------
+
+The required packages are listed in ``environment.yml`` and in ``pyproject.toml``.
+A secondary environment available in  ``dev-environment.yml`` can be used for development, including testing capabilities and jupyter notebooks. 
+
+.. note::
+	Both Unix and MacOS are supported. Python >=3.8 is requested.
+
+
+
+
diff --git a/docs/sphinx/source/intro.rst b/docs/sphinx/source/intro.rst
new file mode 100644
index 0000000..a34eb53
--- /dev/null
+++ b/docs/sphinx/source/intro.rst
@@ -0,0 +1,33 @@
+Introduction
+============
+
+A compact regridder using sparse matrix multiplication
+------------------------------------------------------
+
+`smmregrid` is a tool meant to work into python environment to perform efficiently regridding operation of dask-enable xarray.Dataset and xarray.DataArray.
+It uses sparse matrix multiplication with basic manipulation of the coordinates and uses CDO as a backend to provide the computation of weights. 
+It can work from source and target grids, but most importantly it can also start from precomputed weights.
+
+Please note that `smmregrid` is not meant to be "another interpolation tool", but rather a method to apply pre-computed weights within the python environment. 
+The speedup is estimated to be about ~1.5 to ~5 times compared to CDO itself, slightly lower if then files are written to the disk. 
+
+2D and 3D data are supported on all the grids supported by CDO (including gaussian reduced, healpix or unstructured grids with cell corners), both xarray.Dataset and xarray.DataArray can be used. 
+It relies on the `GridType()` class which identify the properties of each xarray.DataArray, detecting time, horizontal and vertical coordinates.
+
+Masks are treated in a simple way but correctly transfered from source to target grid. 
+
+A specific treatment for mask changing along the vertical dimension can be applied through the `vertical_dim` keyword. 
+Indeed, 3D weights are computed specifically for each level and then stored together in specific files so that it guarantees precise mask handling.
+
+
+.. warning ::
+
+   1. It does not work (yet) with dataset/files including multiple horizontal grids
+   2. It does not work correctly if the Xarray.Dataset includes fields with time-varying missing points
+
+Acknowledgement
+---------------
+
+This repository represents an extension and development of the regridding routines in ``climtas`` by Scott Wales, 
+which already implements efficiently this idea and has no other significant dependencies. 
+
diff --git a/docs/sphinx/source/usage.rst b/docs/sphinx/source/usage.rst
new file mode 100644
index 0000000..7f2f20e
--- /dev/null
+++ b/docs/sphinx/source/usage.rst
@@ -0,0 +1,56 @@
+Basic Usage
+===========
+
+`smmregrid` can be used in two ways. By feeding to the tool a previously generated weights computed
+with cdo or by creating the weights starting from the source and the target grid.
+
+For further details on the different usages, please consult the `demo notebook <https://github.com/jhardenberg/smmregrid/blob/main/demo.ipynb>`_
+
+Starting from weights
+---------------------
+
+You can simply generate the weights with `CdoGenerate()` class and its method `weights()` (or load precomputed weights on your own).
+This will create the weight with the `gen<method>` command from CDO. 
+Then activate the sparse smmregrid `Regridder()` class and run the `regrid()` method on the data you want to regrid.
+
+
+.. code-block:: python 
+
+    from smmregrid import Regridder, CdoGenerate
+    weights = CdoGenerate(filein, target_grid).weights(method='con')
+    interpolator = Regridder(weights=wfield)
+    myfile = interpolator.regrid(xfield[var])
+
+Starting from data
+------------------
+
+You can simply initiate the class with the source file (both a file or a xarray object) and then
+simply operate with the `.regrid()` method
+
+.. code-block:: python
+
+    from smmregrid import Regridder
+    interpolator = Regridder(filein, target_grid)
+    myfile = interpolator.regrid(xfield[var])
+
+
+Area generation
+---------------
+
+Similarly, it is possible to produce the area for both source and target grid making use of 
+the `CdoGenerate()` class and the method `areas()`. This build again on CDO, specifically on
+the `gridarea` command.
+
+For source area:
+
+.. code-block:: python 
+
+    from smmregrid import CdoGenerate
+    areas = CdoGenerate(filein, target_grid).areas()
+
+For target area:
+
+.. code-block:: python 
+
+    from smmregrid import CdoGenerate
+    areas = CdoGenerate(filein, target_grid).areas(target=True)
\ No newline at end of file
diff --git a/environment.yml b/environment.yml
index f5d6593..4b5c901 100644
--- a/environment.yml
+++ b/environment.yml
@@ -7,7 +7,7 @@ name : smmregrid
 channels:
   - conda-forge
 dependencies:
-  - python>=3.8
+  - python>=3.8,<3.13
   - ecCodes
   - numba
   - cdo
diff --git a/pyproject.toml b/pyproject.toml
index 2f75a00..69e5241 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,19 +39,19 @@ tests = [
     "pytest"
 ]
 
+docs = [
+    "sphinx<7.0",
+    "sphinx-rtd-theme"
+]
+
 notebooks = [
     "ipykernel",
     "matplotlib",
     "jinja2"
 ]
 
-#docs = [
-#    "sphinx",
-#    "sphinx-rtd-theme"
-#]
-
 all = [
-    "smmregrid[tests,notebooks]"
+    "smmregrid[tests,notebooks,docs]"
 ]
 
 [project.urls]
diff --git a/smmregrid/__init__.py b/smmregrid/__init__.py
index 72b9d1c..7958324 100644
--- a/smmregrid/__init__.py
+++ b/smmregrid/__init__.py
@@ -1,7 +1,12 @@
 """smmregrid module"""
 
 from .regrid import Regridder, regrid
-from .cdo_weights import cdo_generate_weights
-# from .util import find_vert_coords
+from .cdogenerate import cdo_generate_weights, CdoGenerate
+from .gridtype import GridType, DEFAULT_DIMS
+from .gridinspector import GridInspector
 
-__version__ = '0.0.7'
+
+__version__ = '0.1.0'
+
+__all__ = ['Regridder', 'regrid', 'cdo_generate_weights',
+           'GridType', 'GridInspector', 'CdoGenerate']
diff --git a/smmregrid/cdo_weights.py b/smmregrid/cdo_weights.py
deleted file mode 100644
index e5d4f27..0000000
--- a/smmregrid/cdo_weights.py
+++ /dev/null
@@ -1,271 +0,0 @@
-"""CDO-based generation of weights"""
-
-import os
-import sys
-import tempfile
-import subprocess
-import warnings
-from multiprocessing import Process, Manager
-import numpy
-import xarray
-from .util import find_vert_coords
-from .weights import compute_weights_matrix3d, compute_weights_matrix, mask_weights, check_mask
-from .log import setup_logger
-
-
-
-def worker(wlist, nnn, *args, **kwargs):
-    """Run a worker process"""
-    wlist[nnn] = cdo_generate_weights2d(*args, **kwargs).compute()
-
-
-def cdo_generate_weights(source_grid, target_grid, method="con", extrapolate=True,
-                         remap_norm="fracarea", remap_area_min=0.0, icongridpath=None,
-                         gridpath=None, extra=None, cdo_extra=None, cdo_options=None, vert_coord=None,
-                         cdo="cdo", nproc=1, loglevel='warning'):
-    """Generate the weights using CDO, handling both 2D and 3D cases"""
-
-    loggy = setup_logger(level=loglevel, name='smmregrid.cdo_generate_weights')
-
-    # Check for deprecated 'extra' argument
-    if extra is not None:
-        warnings.warn(
-            "'extra' is deprecated and will be removed in future versions. "
-            "Please use 'cdo_extra' instead.",
-            DeprecationWarning
-        )
-        # If cdo_extra is not provided, use the value from extra
-        if cdo_extra is None:
-            cdo_extra = extra
-
-    # Check if there is a vertical coordinate for 3d oceanic data
-    if not vert_coord:
-        vert_coord = find_vert_coords(source_grid)
-        loggy.info('vert_coord is %s', str(vert_coord))
-
-    if not vert_coord:  # Are we 2D? Use default method
-        weights = cdo_generate_weights2d(
-            source_grid,
-            target_grid,
-            method=method,
-            extrapolate=extrapolate,
-            remap_norm=remap_norm,
-            remap_area_min=remap_area_min,
-            icongridpath=icongridpath,
-            gridpath=gridpath,
-            cdo_extra=cdo_extra,
-            cdo_options=cdo_options,
-            cdo=cdo,
-            nproc=nproc)
-
-        # Precompute destination weights mask
-        weights_matrix = compute_weights_matrix(weights)
-        weights = mask_weights(weights, weights_matrix, vert_coord)
-        masked = int(check_mask(weights, vert_coord))
-        masked_xa = xarray.DataArray(masked, name="dst_grid_masked")
-
-        return xarray.merge([weights, masked_xa])
-
-    else:  # we are 3D
-        cdo_extra = cdo_extra if isinstance(cdo_extra, list) else ([cdo_extra] if cdo_extra else [])
-
-        if isinstance(source_grid, str):
-            sgrid = xarray.open_dataset(source_grid)
-        else:
-            sgrid = source_grid
-
-        nvert = sgrid[vert_coord].values.size
-        # print(nvert)
-
-        # for lev in range(0, nvert):
-        mgr = Manager()
-
-        # dictionaries are shared, so they have to be passed as functions
-        wlist = mgr.list(range(nvert))
-
-        num_blocks, remainder = divmod(nvert, nproc)
-        num_blocks = num_blocks + (0 if remainder == 0 else 1)
-
-        blocks = numpy.array_split(numpy.arange(nvert), num_blocks)
-        for block in blocks:
-            processes = []
-            for lev in block:
-                loggy.info("Generating level: %s", str(lev))
-                cdo_extra_vertical = [f"-sellevidx,{lev+1}"]
-                ppp = Process(target=worker,
-                              args=(wlist, lev, source_grid, target_grid),
-                              kwargs={
-                                    "method": method,
-                                    "extrapolate": extrapolate,
-                                    "remap_norm": remap_norm,
-                                    "remap_area_min": remap_area_min,
-                                    "icongridpath": icongridpath,
-                                    "gridpath": gridpath,
-                                    "cdo_extra": cdo_extra + cdo_extra_vertical,
-                                    "cdo_options": cdo_options,
-                                    "cdo": cdo,
-                                    "nproc": nproc
-                                })
-                ppp.start()
-                processes.append(ppp)
-
-            for proc in processes:
-                proc.join()
-
-        weights = weightslist_to_3d(wlist, vert_coord)
-
-        # Precompute destination weights mask
-        weights_matrix = compute_weights_matrix3d(weights, vert_coord)
-        weights = mask_weights(weights, weights_matrix, vert_coord)
-        masked = check_mask(weights, vert_coord)
-        masked = [int(x) for x in masked]  # convert to list of int
-        masked_xa = xarray.DataArray(masked, 
-                                     coords={vert_coord: range(0, len(masked))},
-                                     name="dst_grid_masked")
-
-        return xarray.merge([weights, masked_xa])
-
-
-def cdo_generate_weights2d(source_grid, target_grid, method="con", extrapolate=True,
-                           remap_norm="fracarea", remap_area_min=0.0, icongridpath=None,
-                           gridpath=None, cdo_extra=None, cdo_options=None, cdo="cdo",
-                           nproc=1):
-    """
-    Generate weights for regridding using CDO
-
-    Available weight generation methods are:
-
-     * bic: SCRIP Bicubic
-     * bil: SCRIP Bilinear
-     * con: SCRIP First-order conservative
-     * con2: SCRIP Second-order conservative
-     * dis: SCRIP Distance-weighted average
-     * laf: YAC Largest area fraction
-     * ycon: YAC First-order conservative
-     * nn: Nearest neighbour
-
-    Run ``cdo gen${method} --help`` for details of each method
-
-    Args:
-        source_grid (xarray.DataArray): Source grid
-        target_grid (xarray.DataArray): Target grid
-            description
-        method (str): Regridding method - default conservative
-        extrapolate (bool): Extrapolate output field
-        remap_norm (str): Normalisation method for conservative methods
-        remap_area_min (float): Minimum destination area fraction
-        gridpath (str): where to store downloaded grids
-        icongridpath (str): location of ICON grids (e.g. /pool/data/ICON)
-        cdo_extra: command(s) to apply to source grid before weight generation (can be a list)
-        cdo_options: command(s) to apply to cdo (can be a list)
-        cdo: the command to launch cdo ["cdo"]
-        nproc: number of processes to use for weight generation (NOT USED!)
-
-        cdo: path to cdo binary
-    Returns:
-        :obj:`xarray.Dataset` with regridding weights
-    """
-
-    supported_methods = ["bic", "bil", "con", "con2", "dis", "laf", "nn", "ycon"]
-    if method not in supported_methods:
-        raise ValueError('The remap method provided is not supported!')
-    if remap_norm not in ["fracarea", "destarea"]:
-        raise ValueError('The remap normalization provided is not supported!')
-
-    # Make some temporary files that we'll feed to CDO
-    weight_file = tempfile.NamedTemporaryFile()
-
-    if isinstance(source_grid, str):
-        sgrid = source_grid
-    else:
-        source_grid_file = tempfile.NamedTemporaryFile()
-        source_grid.to_netcdf(source_grid_file.name)
-        sgrid = source_grid_file.name
-
-    if isinstance(target_grid, str):
-        tgrid = target_grid
-    else:
-        target_grid_file = tempfile.NamedTemporaryFile()
-        target_grid.to_netcdf(target_grid_file.name)
-        tgrid = target_grid_file.name
-
-    # Setup environment
-    env = os.environ
-    if extrapolate:
-        env["REMAP_EXTRAPOLATE"] = "on"
-    else:
-        env["REMAP_EXTRAPOLATE"] = "off"
-
-    env["CDO_REMAP_NORM"] = remap_norm
-    env["REMAP_AREA_MIN"] = f"{remap_area_min:f}"
-
-    if gridpath:
-        env["CDO_DOWNLOAD_PATH"] = gridpath
-    if icongridpath:
-        env["CDO_ICON_GRIDS"] = icongridpath
-
-    try:
-
-        # condense cdo options
-        cdo_extra = cdo_extra if isinstance(cdo_extra, list) else ([cdo_extra] if cdo_extra else [])
-        cdo_options = cdo_options if isinstance(cdo_options, list) else ([cdo_options] if cdo_options else [])
-
-        command = [
-            cdo,
-            *cdo_options,
-            f"gen{method},{tgrid}",  # Method and target grid
-            *cdo_extra,
-            sgrid,
-            weight_file.name
-        ]
-
-        # call to subprocess in compact way
-        subprocess.check_output(
-            command,
-            stderr=subprocess.STDOUT,
-            env=env,
-        )
-
-        # Grab the weights file it outputs as a xarray.Dataset
-        weights = xarray.open_dataset(weight_file.name, engine="netcdf4")
-        return weights
-
-    except subprocess.CalledProcessError as err:
-
-        print(err.output.decode(), file=sys.stderr)
-        raise
-
-    finally:
-        # Clean up the temporary files
-        if not isinstance(source_grid, str):
-            source_grid_file.close()
-        if not isinstance(target_grid, str):
-            target_grid_file.close()
-        weight_file.close()
-
-
-def weightslist_to_3d(ds_list, vert_coord='lev'):
-    """
-    Function to combine a list of 2D cdo weights into a 3D one adding a vertical coordinate lev
-    """
-    # CDO 2.2.0 fix
-    if "numLinks" in ds_list[0].dims:
-        links_dim = "numLinks"
-    else:
-        links_dim = "num_links"
-
-    dim_values = range(len(ds_list))
-    nl = [ds.src_address.size for ds in ds_list]
-    nl0 = max(nl)
-    nlda = xarray.DataArray(nl, coords={vert_coord: range(0, len(nl))}, name="link_length")
-    new_array = []
-    varlist = ["src_address", "dst_address", "remap_matrix", "src_grid_imask", "dst_grid_imask"]
-    ds0 = ds_list[0].drop_vars(varlist)
-    for x, d in zip(ds_list, dim_values):
-        nl1 = x.src_address.size
-#        xplist = [x[vname].pad(num_links=(0, nl0-nl1), mode='constant', constant_values=0)
-        xplist = [x[vname].pad(**{links_dim: (0, nl0 - nl1), "mode": 'constant', "constant_values": 0})
-                  for vname in varlist]
-        xmerged = xarray.merge(xplist)
-        new_array.append(xmerged.assign_coords({vert_coord: d}))
-    return xarray.merge([nlda, ds0, xarray.concat(new_array, vert_coord)])
diff --git a/smmregrid/cdogenerate.py b/smmregrid/cdogenerate.py
new file mode 100644
index 0000000..b8e6136
--- /dev/null
+++ b/smmregrid/cdogenerate.py
@@ -0,0 +1,415 @@
+"""CDO-based generation of weights"""
+
+import os
+import sys
+import copy
+import tempfile
+import subprocess
+import warnings
+from multiprocessing import Process, Manager
+import numpy
+import xarray
+from .weights import compute_weights_matrix3d, compute_weights_matrix, mask_weights, check_mask
+from .log import setup_logger
+from .util import check_gridfile, deprecated_argument
+
+
+class CdoGenerate():
+    """CDO-based class to generate weights for smmregrid."""
+
+    def __init__(self, source_grid, target_grid=None, extra=None, cdo_extra=None,
+            cdo_options=None, cdo_download_path=None, cdo_icon_grids=None,
+            cdo="cdo", loglevel='warning'):
+        """
+        Initialize GenerateWeights class for regridding using Climate Data Operators (CDO), 
+        accommodating both 2D and 3D grid cases.
+
+        Args:
+            source_grid (str or xarray.Dataset): The source grid from which to generate weights.
+                                                 This can be a file path or an xarray dataset.
+            target_grid (str or xarray.Dataset): The target grid to which the source grid 
+                                                 will be regridded. This can also be a file 
+                                                 path or an xarray dataset.
+            loglevel (str, optional): The logging level for messages. Default is 'warning'. 
+                                      Options include 'debug', 'info', 'warning', 
+                                      'error', and 'critical'.
+            extra (list, optional): Deprecated. Previously used for additional CDO options. 
+                                   Use `cdo_extra` instead.
+            cdo_extra (list, optional): Additional CDO command-line options. 
+                                               Defaults to None.
+            cdo_options (list, optional): Options for CDO commands. Defaults to None.
+            cdo (str, optional): The command to invoke CDO. Default is "cdo".
+            cdo_icon_grids (str, optional): Path to the ICON grid 
+                                            if applicable. Defaults to None.
+            cdo_download_path (str, optional): Path to the grid download path 
+                                                if applicable. Defaults to None.
+        """
+
+        self.loggy = setup_logger(level=loglevel, name='smmregrid.CdoGenerate')
+
+        # deprecation warning
+        cdo_extra = deprecated_argument(extra, cdo_extra, 'extra', 'cdo_extra')
+
+        # cdo options and extra, ensure they are lists
+        self.cdo_extra = cdo_extra if isinstance(cdo_extra, list) else ([cdo_extra] if cdo_extra else [])
+        self.cdo_options = cdo_options if isinstance(cdo_options, list) else ([cdo_options] if cdo_options else [])
+
+        # assign the two grids
+        self.source_grid = source_grid
+        self.target_grid = target_grid
+
+        # define grid filenames
+        self.source_grid_filename = None
+        self.target_grid_filename = None
+
+        # assign cdo bin file and get envrinment file
+        self.cdo = cdo
+        self.env = os.environ.copy()
+        if cdo_download_path:
+            self.env["CDO_DOWNLOAD_PATH"] = cdo_download_path
+        if cdo_icon_grids:
+            self.env["CDO_ICON_GRIDS"] = cdo_icon_grids
+
+    @staticmethod
+    def _safe_check(method, remap_norm):
+        """Safety checks for weights generation """
+
+        supported_methods = ["bic", "bil", "con", "con2", "dis", "laf", "nn", "ycon"]
+        if method not in supported_methods:
+            raise ValueError('The remap method provided is not supported!')
+        if remap_norm not in ["fracarea", "destarea"]:
+            raise ValueError('The remap normalization provided is not supported!')
+    
+    @staticmethod
+    def _prepare_grid(grid):
+        """Helper function to prepare grid (file or dataset)."""
+
+        if isinstance(grid, str):
+            return grid
+
+        grid_file = tempfile.NamedTemporaryFile(delete=False)
+        grid.to_netcdf(grid_file.name)
+        return grid_file.name
+
+    def weights(self, method="con", extrapolate=True,
+            remap_norm="fracarea", remap_area_min=0.0,
+            vert_coord=None, vertical_dim=None, nproc=1):
+        """
+        Generate weights for regridding using Climate Data Operators (CDO), 
+        accommodating both 2D and 3D grid cases.
+
+        Args:
+            method (str, optional): The remapping method to use. 
+                                    Default is "con" for conservative remapping.
+                                    Other options may include 'bil', 'nearest', etc.
+            extrapolate (bool, optional): Whether to allow extrapolation beyond the grid boundaries. 
+                                          Defaults to True.
+            remap_norm (str, optional): The normalization method to apply when remapping.
+                                        Default is "fracarea" which normalizes by fractional area.
+            remap_area_min (float, optional): Minimum area for remapping. Defaults to 0.0.
+            nproc (int, optional): Number of processes to use for parallel processing. Default is 1.
+            vertical_dim (str, optional): Name of the vertical dimension in the source grid, if applicable.
+                                        Defaults to None. Use if the grid is 3D.
+            vert_coord (str, optional): Deprecated. Previously used to specify the vertical coordinate.
+                                        Use `vertical_dim` instead.
+
+        Returns:
+            xarray.Dataset: A dataset containing the generated weights 
+                            and a mask indicating which grid cells
+                            were successfully masked. 
+                            The mask is stored in a variable named "dst_grid_masked".
+
+        Raises:
+            KeyError: If the specified vertical dimension cannot be found in the source grid.
+            Warning: If deprecated arguments `extra` or `vert_coord` are used.
+
+        Notes:
+            This function handles both 2D and 3D grid cases:
+
+            For 2D grids (when `vertical_dim` is None), it calls the `_weights_2d` method
+            to generate weights. The weights are then masked based on a precomputed weights matrix.
+
+            For 3D grids (when `vertical_dim` is specified), it uses multiprocessing to generate weights
+            for each vertical level. It requires the vertical dimension to be present in the source grid,
+            and it will generate a mask indicating valid and invalid weights for each vertical level.
+
+            The function logs the progress of weight generation, including the length of vertical dimensions
+            and each level being processed.
+
+            Deprecation Warning: The `extra` and `vert_coord` parameters are deprecated and will be removed in future versions.
+            Users should migrate to using `cdo_extra` and `vertical_dim`, respectively.
+
+        """
+        if self.target_grid is None:
+            raise TypeError('Target grid is not specified, cannot provide any regridding')
+        if self.source_grid is None:
+            raise TypeError('Source grid is not specified, cannot provide any regridding')
+
+        # verify that method and normalization are suitable
+        self._safe_check(method, remap_norm)
+
+        # vertical dimension
+        vertical_dim = deprecated_argument(vert_coord, vertical_dim, 'vert_coord', 'vertical_dim')
+
+        # set grids filenames
+        self.source_grid_filename = self._prepare_grid(self.source_grid)
+        self.target_grid_filename = self._prepare_grid(self.target_grid)
+
+        # Generate weights for 2D or 3D grid based on vertical_dim presence
+        if not vertical_dim:
+            return self._weights_2d(method, extrapolate, remap_norm, remap_area_min)
+        return self._weights_3d(method, extrapolate, remap_norm, remap_area_min,
+                            nproc, vertical_dim)
+
+    def _weights_2d(self, method, extrapolate, remap_norm,
+                remap_area_min):
+        """Generate 2D weights using CDO."""
+
+        weights = self._cdo_generate_weights(method, extrapolate,
+                                             remap_norm, remap_area_min)
+        weights_matrix = compute_weights_matrix(weights)
+        weights = mask_weights(weights, weights_matrix)
+        masked = int(check_mask(weights))
+        masked_xa = xarray.DataArray(masked, name="dst_grid_masked")
+        return xarray.merge([weights, masked_xa])
+
+    def _weights_3d(self, method, extrapolate, remap_norm,
+                remap_area_min, nproc, vertical_dim):
+        """Generate 3D weights using multiprocessing."""
+
+        if isinstance(self.source_grid, str):
+            sgrid = xarray.open_dataset(self.source_grid)
+        else:
+            sgrid = self.source_grid
+
+        if vertical_dim not in sgrid.dims:
+            raise KeyError(f'Cannot find vertical dim {vertical_dim} in {list(sgrid.dims)}')
+
+        #nvert = sgrid[vertical_dim].values.size
+        nvert = sgrid.sizes[vertical_dim]
+        self.loggy.info('Vertical dimension has length: %s', nvert)
+
+        mgr = Manager()
+        wlist = mgr.list(range(nvert))
+
+        num_blocks, remainder = divmod(nvert, nproc)
+        num_blocks += 0 if remainder == 0 else 1
+
+        blocks = numpy.array_split(numpy.arange(nvert), num_blocks)
+        for block in blocks:
+            processes = []
+            for lev in block:
+                self.loggy.info("Generating level: %s", str(lev))
+                cdo_extra_vertical = [f"-sellevidx,{lev + 1}"]
+                ppp = Process(target=self._weights_worker,
+                              args=(wlist, lev),
+                              kwargs={
+                                  "method": method,
+                                  "extrapolate": extrapolate,
+                                  "remap_norm": remap_norm,
+                                  "remap_area_min": remap_area_min,
+                                  "cdo_extra_vertical": cdo_extra_vertical
+                              })
+                ppp.start()
+                processes.append(ppp)
+
+            for proc in processes:
+                proc.join()
+
+        weights = self.weightslist_to_3d(wlist, vertical_dim)
+        weights_matrix = compute_weights_matrix3d(weights, vertical_dim)
+        weights = mask_weights(weights, weights_matrix, vertical_dim)
+        masked = check_mask(weights, vertical_dim)
+        masked = [int(x) for x in masked]
+        masked_xa = xarray.DataArray(masked,
+                                     coords={vertical_dim: range(0, len(masked))},
+                                     name="dst_grid_masked")
+
+        return xarray.merge([weights, masked_xa])
+
+    def _weights_worker(self, wlist, nnn, *args, **kwargs):
+        """Run a worker process."""
+        wlist[nnn] = self._cdo_generate_weights(*args, **kwargs).compute()
+
+    def _cdo_generate_weights(self, method="con", extrapolate=True,
+                              remap_norm="fracarea", remap_area_min=0.0,
+                              cdo_extra_vertical=None):
+        """
+        Generate weights for regridding using CDO
+
+        Available weight generation methods are:
+
+        * bic: SCRIP Bicubic
+        * bil: SCRIP Bilinear
+        * con: SCRIP First-order conservative
+        * con2: SCRIP Second-order conservative
+        * dis: SCRIP Distance-weighted average
+        * laf: YAC Largest area fraction
+        * ycon: YAC First-order conservative
+        * nn: Nearest neighbour
+
+        Args:
+            method (str): Regridding method - default "con"
+            extrapolate (bool): Extrapolate output field - default True
+            remap_norm (str): Normalisation method for conservative methods
+            remap_area_min (float): Minimum destination area fraction
+            cdo_extra_vertical (list): Command to select vertical levels for 3d case
+
+        Returns:
+            :obj:`xarray.Dataset` with regridding weights
+        """
+
+        sgrid = self.source_grid_filename
+        tgrid = self.target_grid_filename
+
+        cdo_extra_vertical = cdo_extra_vertical if isinstance(cdo_extra_vertical, list) else ([cdo_extra_vertical] if cdo_extra_vertical else [])
+
+        # Log method and remapping information
+        self.loggy.info("CDO remapping method: %s", method)
+        self.loggy.info("Extrapolation enabled: %s", extrapolate)
+        self.loggy.debug("Normalization method: %s", remap_norm)
+        self.loggy.debug("Minimum remap area: %s", remap_area_min)
+
+        weight_file = tempfile.NamedTemporaryFile()
+        self.loggy.debug("Weight file name is: %s", weight_file.name)
+
+        self.loggy.debug("Source grid file name is: %s", sgrid)
+        self.loggy.debug("Target grid file name is: %s", tgrid)
+
+        env = copy.deepcopy(self.env)
+        env["REMAP_EXTRAPOLATE"] = "on" if extrapolate else "off"
+        env["CDO_REMAP_NORM"] = remap_norm
+        env["REMAP_AREA_MIN"] = f"{remap_area_min:f}"
+
+        try:
+            self.loggy.info("Additional CDO commands: %s", self.cdo_extra)
+            self.loggy.info("Additional CDO options: %s", self.cdo_options)
+
+            command = [
+                self.cdo,
+                *self.cdo_options,
+                f"gen{method},{tgrid}",
+                *self.cdo_extra+cdo_extra_vertical,
+                sgrid,
+                weight_file.name
+            ]
+            self.loggy.debug("Final CDO command: %s", command)
+            subprocess.check_output(command, stderr=subprocess.STDOUT, env=env)
+
+            weights = xarray.open_dataset(weight_file.name, engine="netcdf4")
+            return weights
+    
+        except subprocess.CalledProcessError as err:
+            print(err.output.decode(), file=sys.stderr)
+            raise
+        finally:
+            weight_file.close()
+
+    # def _remove_tmpfile(self, tmpfile):
+    #     """Helper function to clean the tempfile grid (file or dataset)."""
+    #     if not isinstance(tmpfile, str):
+    #         os.remove(tmpfile.name)
+
+    def weightslist_to_3d(self, ds_list, vertical_dim='lev'):
+        """Combine a list of 2D CDO weights into a 3D one."""
+
+        links_dim = "numLinks" if "numLinks" in ds_list[0].dims else "num_links"
+        dim_values = range(len(ds_list))
+        nl = [ds.src_address.size for ds in ds_list]
+        nl0 = max(nl)
+        nlda = xarray.DataArray(nl, coords={vertical_dim: range(0, len(nl))}, name="link_length")
+
+        new_array = []
+        varlist = ["src_address", "dst_address", "remap_matrix", "src_grid_imask", "dst_grid_imask"]
+        ds0 = ds_list[0].drop_vars(varlist)
+
+        for x, d in zip(ds_list, dim_values):
+            nl1 = x.src_address.size
+            xplist = [x[vname].pad(**{links_dim: (0, nl0 - nl1), "mode": 'constant', "constant_values": 0})
+                      for vname in varlist]
+            xmerged = xarray.merge(xplist)
+            new_array.append(xmerged.assign_coords({vertical_dim: d}))
+
+        return xarray.merge([nlda, ds0, xarray.concat(new_array, vertical_dim)],
+                            combine_attrs='no_conflicts')
+    
+    def areas(self, target=False):
+        """Generate source areas or target areas"""
+
+        if not target:
+            self.loggy.info('Generating areas for source grid!')
+            #if not self.source_grid:
+            #    raise TypeError('Source grid is not specified, cannot provide any area')
+            return self._areas(self.source_grid, cdo_extra=self.cdo_extra,
+                                cdo_options=self.cdo_options)
+
+        if self.target_grid:
+            
+            self.loggy.info('Generating areas for target grid!')
+            return self._areas(self.target_grid)
+
+        raise TypeError('Target grid is not specified, cannot provide any area')
+
+            
+    def _areas(self, filename, cdo_extra=None, cdo_options=None):
+        """Generate areas in a similar way of what done for weights"""
+
+        # safety listing
+        cdo_extra = cdo_extra if isinstance(cdo_extra, list) else ([cdo_extra] if cdo_extra else [])
+        cdo_options = cdo_options if isinstance(cdo_options, list) else ([cdo_options] if cdo_options else [])
+
+        # Make some temporary files that we'll feed to CDO
+        areas_file = tempfile.NamedTemporaryFile()
+
+        # prepare grid
+        if check_gridfile(filename) == 'grid':
+            self.loggy.info('CDO grid as %s to be used for area generation', filename)
+            sgrid = f"-const,1,{filename}"
+        else:
+            sgrid = self._prepare_grid(filename)
+
+        try:
+            self.loggy.info("Additional CDO commands: %s", cdo_extra)
+            self.loggy.info("Additional CDO options: %s", cdo_options)
+
+            command = [
+                self.cdo,
+                *cdo_options + ["-f", "nc4"],
+                "gridarea",
+                *cdo_extra,
+                sgrid,
+                areas_file.name
+            ]
+            self.loggy.debug("Final CDO command: %s", command)
+            subprocess.check_output(command, stderr=subprocess.STDOUT, env=self.env)
+
+            areas = xarray.open_dataset(areas_file.name, engine="netcdf4")
+            areas.cell_area.attrs['units'] = 'm2'
+            areas.cell_area.attrs['standard_name'] = 'area'
+            areas.cell_area.attrs['long_name'] = 'area of grid cell'
+            return areas
+    
+        except subprocess.CalledProcessError as err:
+            print(err.output.decode(), file=sys.stderr)
+            raise
+        finally:
+            areas_file.close()
+
+
+def cdo_generate_weights(source_grid, target_grid, method="con", extrapolate=True,
+                         remap_norm="fracarea", remap_area_min=0.0, gridpath=None,
+                         icongridpath=None, extra=None, cdo_extra=None, cdo_options=None,
+                         vertical_dim=None, vert_coord=None,
+                         cdo="cdo", nproc=1, loglevel='warning'):
+    """
+    Wrapper function for the new cdo class to provide the old access
+    """
+    warnings.warn("cdo_generate_weights() is now deprecated, please use CdoGenerate().weights()",
+                    DeprecationWarning)
+    generator = CdoGenerate(source_grid=source_grid, target_grid=target_grid, loglevel=loglevel,
+                            extra=extra, cdo_extra=cdo_extra, cdo_options=cdo_options, cdo=cdo, 
+                            cdo_icon_grids=icongridpath, cdo_download_path=gridpath)
+    return generator.weights(method=method, extrapolate=extrapolate,
+                         remap_norm=remap_norm, remap_area_min=remap_area_min,
+                         vertical_dim=vertical_dim, vert_coord=vert_coord, nproc=nproc)
+
diff --git a/smmregrid/checker.py b/smmregrid/checker.py
index 289eed4..32e98c9 100644
--- a/smmregrid/checker.py
+++ b/smmregrid/checker.py
@@ -3,9 +3,7 @@
 import numpy as np
 import xarray as xr
 from cdo import Cdo
-from smmregrid import Regridder, cdo_generate_weights
-
-
+from smmregrid import Regridder, CdoGenerate
 cdo = Cdo()
 
 
@@ -24,30 +22,38 @@ def find_var(xfield):
 
 
 def check_cdo_regrid(finput, ftarget, remap_method='con', access='Dataset',
-                     init_method='grids', vert_coord=None):
+                     init_method='grids', vertical_dim=None, extrapolate=True,
+                     loglevel='INFO'):
     """Given a file to be interpolated finput over the ftarget grid,
     check if the output of the last variable is the same as produced
     by CDO remap command. This function is used for tests."""
 
     # define files and open input file
-    xfield = xr.open_mfdataset(finput)
+    if isinstance(finput, str):
+        xfield = xr.open_mfdataset(finput)
+    else:
+        xfield = finput
+
+    # convert extrapolation to cdo!
+    cdoextrapolate = 'on' if extrapolate else 'off'
 
     # var as the last available
     # myvar = list(xfield.data_vars)[-1]
 
     # interpolation with pure CDO
     cdo_interpolator = getattr(cdo, 'remap' + remap_method)
-    cdofield = cdo_interpolator(ftarget, input=finput, returnXDataset=True)
+    cdofield = cdo_interpolator(ftarget, input=finput, returnXDataset=True,
+                                env={'REMAP_EXTRAPOLATE': cdoextrapolate})
     # print(cdofield)
 
     # var as the one which have time and not have bnds (could work)
     smmvar = find_var(xfield)
     cdovar = find_var(cdofield)
 
-    if len(smmvar) == 1 and access == 'DataArray':
-        xfield = xfield[smmvar[0]]
-    if len(cdovar) == 1 and access == 'DataArray':
-        cdofield = cdofield[cdovar[0]]
+    # if len(smmvar) == 1 and access == 'DataArray':
+    #    xfield = xfield[smmvar[0]]
+    # if len(cdovar) == 1 and access == 'DataArray':
+    #    cdofield = cdofield[cdovar[0]]
 
     # interpolation with smmregrid (CDO-based)
     # method with creation of weights
@@ -57,11 +63,12 @@ def check_cdo_regrid(finput, ftarget, remap_method='con', access='Dataset',
     # method with automatic creation of weights
     if init_method == 'grids':
         interpolator = Regridder(source_grid=finput, target_grid=ftarget,
-                                 method=remap_method, vert_coord=vert_coord)
-    if init_method == 'weights':
-        wfield = cdo_generate_weights(finput, ftarget,
-                                      method=remap_method, vert_coord=vert_coord)
-        interpolator = Regridder(weights=wfield, vert_coord=vert_coord)
+                                 method=remap_method, vertical_dim=vertical_dim, loglevel=loglevel)
+    elif init_method == 'weights':
+        wfield = CdoGenerate(finput, ftarget, loglevel=loglevel).weights(method=remap_method, vertical_dim=vertical_dim)
+        interpolator = Regridder(weights=wfield, loglevel=loglevel)
+    else:
+        raise KeyError('Unsupported init method')
     rfield = interpolator.regrid(xfield)
 
     if access == 'Dataset':
@@ -73,7 +80,8 @@ def check_cdo_regrid(finput, ftarget, remap_method='con', access='Dataset',
     return checker
 
 
-def check_cdo_regrid_levels(finput, ftarget, vert_coord, levels, remap_method='con', access='Dataset'):
+def check_cdo_regrid_levels(finput, ftarget, vertical_dim, levels, remap_method='con',
+                            access='Dataset', extrapolate=True, loglevel='INFO'):
     """Given a file to be interpolated finput over the ftarget grid,
     check if the output of the last variable is the same as produced
     by CDO remap command. This function is used for tests.
@@ -82,35 +90,30 @@ def check_cdo_regrid_levels(finput, ftarget, vert_coord, levels, remap_method='c
     # define files and open input file
     xfield = xr.open_mfdataset(finput)
 
+    # convert extrapolate to CDO
+    cdoextrapolate = 'on' if extrapolate else 'off'
+
     # interpolation with pure CDO
     cdo_interpolator = getattr(cdo, 'remap' + remap_method)
-    cdofield = cdo_interpolator(ftarget, input=finput, returnXDataset=True)
+    cdofield = cdo_interpolator(ftarget, input=finput, returnXDataset=True,
+                                env={'REMAP_EXTRAPOLATE': cdoextrapolate})
 
     # Keep only some levels
-    cdofield = cdofield.isel(**{vert_coord: levels})
+    cdofield = cdofield.isel(**{vertical_dim: levels})
 
     # var as the one which have time and not have bnds (could work)
     smmvar = find_var(xfield)
     cdovar = find_var(cdofield)
 
-    if len(smmvar) == 1 and access == 'DataArray':
-        xfield = xfield[smmvar[0]]
-    if len(cdovar) == 1 and access == 'DataArray':
-        cdofield = cdofield[cdovar[0]]
-
     # compute weights
-    wfield = cdo_generate_weights(finput, ftarget,
-                                  method=remap_method, vert_coord=vert_coord)
-    
-    # Pass full 3D weights
-    interpolator = Regridder(weights=wfield, vert_coord=vert_coord)
+    wfield = CdoGenerate(finput, ftarget, loglevel=loglevel).weights(
+        method=remap_method, vertical_dim=vertical_dim)
 
-    # Add a helper idx_3d coordinate
-    idx = list(range(0, len(xfield.coords[vert_coord])))
-    xfield = xfield.assign_coords(idx_3d=(vert_coord, idx))
+    # Pass full 3D weights
+    interpolator = Regridder(weights=wfield, loglevel=loglevel)
 
     # subselect some levels
-    xfield = xfield.isel(**{vert_coord: levels})
+    xfield = xfield.isel(**{vertical_dim: levels})
 
     # Regrid level selection
     rfield = interpolator.regrid(xfield)
@@ -121,4 +124,4 @@ def check_cdo_regrid_levels(finput, ftarget, vert_coord, levels, remap_method='c
 
     # check if arrays are equal with numerical tolerance
     checker = np.allclose(cdofield, rfield, equal_nan=True)
-    return checker
\ No newline at end of file
+    return checker
diff --git a/smmregrid/gridinspector.py b/smmregrid/gridinspector.py
new file mode 100644
index 0000000..5ad9a7a
--- /dev/null
+++ b/smmregrid/gridinspector.py
@@ -0,0 +1,194 @@
+"""GridInspector class module"""
+
+import xarray as xr
+from smmregrid.log import setup_logger
+from .gridtype import GridType
+
+
+class GridInspector():
+    """Class to investigate data and detect its GridType() object"""
+
+    def __init__(self, data, cdo_weights=False, extra_dims=None,
+                 clean=True, loglevel='warning'):
+        """
+        GridInspector class to detect information on the data, based on GridType class
+
+        Parameters:
+            data (xr.Datase or xr.DataArray): The input dataset.
+            clean (bool): apply the cleaning of grids which are assumed to be not relevant 
+                          for regridding purposes (e.g. bounds)
+            cdo_weights (bool): if the data provided are cdo weights instead of data to be regridded
+            extra_dims(dict): Extra dimensions to be added and passed to GridType
+            loglevel: The loglevel that you want you use
+        """
+        self.loggy = setup_logger(name='smmregrid.GridInspect', level=loglevel)
+        self.loglevel = loglevel
+        self.extra_dims = extra_dims
+        self.data = data
+        self.cdo_weights = cdo_weights
+        self.clean = clean
+        self.grids = []  # List to hold all grids info
+        self.loggy.debug('Extra_dims are %s', extra_dims)
+        self.loggy.debug('Clean flag is %s and cdo_weights init is %s', clean, cdo_weights)
+
+    def _inspect_grids(self):
+        """
+        Inspects the dataset and identifies different grids.
+        """
+
+        if isinstance(self.data, xr.Dataset):
+            for variable in self.data.data_vars:
+                self._inspect_dataarray_grid(self.data[variable])
+        elif isinstance(self.data, xr.DataArray):
+            self._inspect_dataarray_grid(self.data)
+        else:
+            raise ValueError('Data supplied is neither xarray Dataset or DataArray')
+
+        # get variables associated to the grid
+        for gridtype in self.grids:
+            self.identify_variables(gridtype)
+
+    def _inspect_dataarray_grid(self, data_array):
+        """
+        Helper method to inspect a single DataArray and identify its grid type.
+        """
+        grid_key = tuple(data_array.dims)
+        gridtype = GridType(dims=grid_key, extra_dims=self.extra_dims)
+        if gridtype not in self.grids:
+            self.grids.append(gridtype)
+
+    def _inspect_weights(self):
+        """
+        Return basic information about CDO weights
+        """
+
+        gridtype = GridType(dims=[], weights=self.data)
+
+        # get vertical info from the weights coords if available
+        if self.data.coords:
+            self.loggy.debug('Vertical dimension read from weights and assigned to %s',
+                             list(self.data.coords)[0])
+            gridtype.vertical_dim = list(self.data.coords)[0]
+
+        self.grids.append(gridtype)
+
+    def get_grid_info(self):
+        """
+        Returns detailed information about all the grids in the dataset.
+        """
+
+        if isinstance(self.data, str):
+            self.data = xr.open_dataset(self.data)
+
+        if self.cdo_weights:
+            self.loggy.info('CDO weights are used to define the grid')
+            self._inspect_weights()
+        else:
+            self._inspect_grids()
+
+        if self.clean:
+            self._clean_grids()
+
+        # self.loggy.info('Grids that have been identifed are: %s', self.grids.)
+        for gridtype in self.grids:
+            self.loggy.debug('More details on gridtype %s:', gridtype.dims)
+            if gridtype.horizontal_dims:
+                self.loggy.debug('    Space dims are: %s', gridtype.horizontal_dims)
+            if gridtype.vertical_dim:
+                self.loggy.debug('    Vertical dims is: %s', gridtype.vertical_dim)
+            self.loggy.debug('    Variables are: %s', list(gridtype.variables.keys()))
+            self.loggy.debug('    Bounds are: %s', gridtype.bounds)
+        return self.grids
+
+    def _clean_grids(self):
+        """
+        Remove degenerate grids which are used by not relevant variables
+        """
+        removed = []
+
+        # Iterate through grids
+        for gridtype in self.grids:
+            # Check if any variable in the grid contains 'bnds'
+            if not gridtype.dims:
+                removed.append(gridtype)
+                self.loggy.info('Removing the grid defined by %s with with no spatial dimensions',
+                                gridtype.dims)
+            elif all('bnds' in variable for variable in gridtype.variables):
+                removed.append(gridtype)  # Add to removed list
+                self.loggy.info('Removing the grid defined by %s with variables containing "bnds"',
+                                gridtype.dims)
+            elif all('bounds' in variable for variable in gridtype.variables):
+                removed.append(gridtype)  # Add to removed list
+                self.loggy.info('Removing the grid defined by %s with variables containing "bounds"',
+                                gridtype.dims)
+
+        for remove in removed:
+            self.grids.remove(remove)
+
+
+    def _identify_variable(self, gridtype, var_data, var_name=None):
+        """Helper function to process individual variables.
+
+        Args:
+            gridtype (GridType): The Gridtype object of originally inspected from the data
+            var_data (xr.DataArray): An xarray DataArray containing variable data.
+            var_name (str, optional): The name of the variable. If None, uses the name from var_data.
+
+        Updates:
+            self.variables: Updates the variables dictionary with the variable's coordinates.
+        """
+        datagridtype = GridType(var_data.dims, extra_dims=self.extra_dims)
+
+        if datagridtype == gridtype:
+            gridtype.variables[var_name or var_data.name] = {
+                'coords': list(var_data.coords),
+            }
+
+    def identify_variables(self, gridtype):
+        """
+        Identify variables in the provided data that match the defined dimensions.
+
+        Args:
+            gridtype (GridType): The Gridtype object of originally inspected from the data
+
+        Raises:
+            TypeError: If the input data is neither an xarray Dataset nor DataArray.
+
+        Updates:
+            self.variables: Updates the variables dictionary with identified variables and their coordinates.
+            self.bounds: Updates the bounds list with identified bounds variables from the dataset.
+        """
+
+        if not isinstance(self.data, (xr.Dataset, xr.DataArray)):
+            raise TypeError("Unsupported data type. Must be an xarray Dataset or DataArray.")
+
+        if isinstance(self.data, xr.Dataset):
+            for var in self.data.data_vars:
+                self._identify_variable(gridtype, self.data[var], var)
+            gridtype.bounds = self._identify_spatial_bounds(self.data)
+            gridtype.variables = {
+                            key: value
+                            for key, value in gridtype.variables.items()
+                            if key not in set(gridtype.bounds)
+                        }
+
+        elif isinstance(self.data, xr.DataArray):
+            self._identify_variable(gridtype, self.data)
+
+    def _identify_spatial_bounds(self, data):
+        """
+        Identify bounds variables in the dataset by checking variable names.
+
+        Args:
+            data (xr.Dataset): An xarray dataset containing data variables.
+
+        Returns:
+            list: A list of bounds variable names identified in the dataset.
+        """
+        bounds_variables = []
+
+        for var in data.data_vars:
+            if (var.endswith('_bnds') or var.endswith('_bounds') or var == 'vertices') and 'time' not in var:
+                bounds_variables.append(var)
+
+        return bounds_variables
diff --git a/smmregrid/gridtype.py b/smmregrid/gridtype.py
new file mode 100644
index 0000000..7df27dc
--- /dev/null
+++ b/smmregrid/gridtype.py
@@ -0,0 +1,154 @@
+# GridType class to gather all information about grids with shared dimensions
+
+# default spatial dimensions and vertical coordinates
+DEFAULT_DIMS = {
+    'horizontal': ['i', 'j', 'x', 'y', 'lon', 'lat', 'longitude', 'latitude',
+                   'cell', 'cells', 'ncells', 'values', 'value', 'nod2', 'pix', 'elem',
+                   'nav_lon', 'nav_lat'],
+    'vertical': ['lev', 'nz1', 'nz', 'depth', 'depth_full', 'depth_half'],
+    'time': ['time']
+}
+
+
+class GridType:
+    """Fundamental GridType object"""
+
+    def __init__(self, dims, extra_dims=None, weights=None):
+        """
+        Initializes a GridType object carrying grid-specific information required by smmregrid.
+
+        Args:
+            dims (list): A list of default dimensions for the grid (e.g., ['time', 'lat', 'lon']).
+            extra_dims (dict, optional): A dictionary including keys 'vertical', 'time', and 'horizontal'
+                                          that can be used to extend the default dimensions. Defaults to None.
+            weights (any, optional): CDO weights used in regridding. It will initiate the object in a different
+                                         way assuming single-gridtype objects. Defaults to None.
+
+        Attributes:
+            horizontal_dims (list): The identified horizontal dimensions from the input.
+            vertical_dim (str or None): The identified vertical dimension, if applicable.
+            dims (list): The dimensions defined for the grid. A combination of horizontal and vertical. 
+            time_dims (list): The identified time dimensions from the input.
+            other_dims (list): The dimensions which are there but are not identified automatically.
+            variables (dict): A dictionary holding identified variables and their coordinates.
+            bounds (list): A list of bounds variables identified in the dataset.
+            masked (any): Placeholder for masked data (to be defined later).
+            weights (any): The weights associated with the grid, if provided.
+            weights_matrix (any): Placeholder for a weights matrix (to be defined later).
+            level_index (str): A string used to identify the index of the levels.
+
+        Raises:
+            ValueError: If multiple vertical dimensions are identified during initialization.
+        """
+
+        # key definitions
+        #if extra_dims:
+        #    for key, value in extra_dims.items():
+        #        if not isinstance(value, list):
+        #            raise ValueError(f'Extra_dims for {key} must be a list!')
+
+        default_dims = self._handle_default_dimensions(extra_dims)
+        self.horizontal_dims = self._identify_dims('horizontal', dims, default_dims)
+        self.vertical_dim = self._identify_dims('vertical', dims, default_dims)
+        self.dims = (self.horizontal_dims or []) + ([self.vertical_dim] if self.vertical_dim else [])
+        self.time_dims = self._identify_dims('time', dims, default_dims)
+        self.other_dims = self._identify_other_dims(dims)
+
+        # used by GridInspector
+        self.variables = {}
+        self.bounds = []
+
+        # used by Regrid class
+        self.masked = None
+        self.weights = weights
+        self.weights_matrix = None
+        self.level_index = "idx_"
+
+    def _handle_default_dimensions(self, extra_dims):
+        """
+        Extend the default dimensions based on the provided extra dimensions.
+
+        Args:
+            extra_dims (dict): A dictionary that can include 'vertical', 'time', and
+                               'horizontal' keys for extending dimensions.
+
+        Returns:
+            dict: An updated dictionary of default dimensions that includes any extensions specified
+                  in extra_dims.
+
+        Notes:
+            If extra_dims is None, the default dimensions remain unchanged.
+        """
+
+        if extra_dims is None:
+            return DEFAULT_DIMS
+
+        update_dims = DEFAULT_DIMS
+        for dim in extra_dims.keys():
+            if extra_dims[dim]:
+                update_dims[dim] = list(set(update_dims[dim] + extra_dims[dim]))
+        return update_dims
+
+    def __eq__(self, other):
+        """
+        Check for equality between two GridType instances.
+
+        Args:
+            other (GridType): Another GridType instance to compare against.
+
+        Returns:
+            bool: True if both self.dims are equal for both instances, False otherwise.
+        """
+        if isinstance(other, GridType):
+            return set(self.dims) == set(other.dims)
+        return False
+
+    def __hash__(self):
+        """
+        Generate a hash value for the GridType instance based on its dimensions.
+
+        Returns:
+            int: A hash value representing the dimensions of the GridType instance.
+        """
+        return hash(tuple(self.dims))
+
+    def _identify_dims(self, axis, dims, default_dims):
+        """
+        Identify dimensions along a specified axis.
+
+        Args:
+            axis (str): The axis to check ('horizontal', 'vertical', or 'time').
+            default_dims (dict): The dictionary of default dimensions to check against.
+
+        Returns:
+            list or str: A list of identified dimensions or a single identified vertical dimension.
+                          Returns None if no dimensions are identified.
+
+        Raises:
+            ValueError: If more than one vertical dimension is identified.
+        """
+        identified_dims = list(set(dims).intersection(default_dims[axis]))
+        if axis == 'vertical':
+            if len(identified_dims) > 1:
+                raise ValueError(f'Only one vertical dimension can be processed at the time: check {identified_dims}')
+            if len(identified_dims) == 1:
+                identified_dims = identified_dims[0]  # unlist the single vertical dimension
+        return identified_dims if identified_dims else None
+
+    def _identify_other_dims(self, dims):
+        """
+        Calculate and return the dimensions that are not part of horizontal_dims, vertical_dim, or time_dims.
+
+        Returns:
+            set: A set of unused dimensions that are not part of horizontal, vertical, or time dimensions.
+        """
+        # Safely convert identified dimensions to sets, handling None values
+        horizontal_dims = set(self.horizontal_dims) if self.horizontal_dims else set()
+        vertical_dim = {self.vertical_dim} if self.vertical_dim else set()  # Use {} for single value
+        time_dims = set(self.time_dims) if self.time_dims else set()
+
+        # Return the unused dimensions by subtracting used dimensions from all dimensions
+        return list(set(dims) - (horizontal_dims | vertical_dim | time_dims))
+    
+
+
diff --git a/smmregrid/log.py b/smmregrid/log.py
index 884ddb5..f2bb55b 100644
--- a/smmregrid/log.py
+++ b/smmregrid/log.py
@@ -13,7 +13,7 @@ def setup_logger(level=None, name=None):
     logger = logging.getLogger(name)  # Create a logger specific to your module
 
     if logger.handlers:
-        logger.warning('Logging is already setup with name %s', name)
+        # logger.warning('Logging is already setup with name %s', name)
         if level != logging.getLevelName(logger.getEffectiveLevel()):
             logger.setLevel(loglev)
             logger.info('Updating the log_level to %s', loglev)
@@ -30,12 +30,13 @@ def setup_logger(level=None, name=None):
 
     # Create a handler for the logger
     handler = logging.StreamHandler()
-    #handler.setLevel(loglev)  # Set the desired log level for the handler
+    # handler.setLevel(loglev)  # Set the desired log level for the handler
     handler.setFormatter(formatter)  # Assign the formatter to the handler
     logger.addHandler(handler)  # Add the handler to the logger
 
     return logger
 
+
 def convert_logger(loglev=None):
     """Convert a string or integer to a valid logging level"""
 
@@ -67,4 +68,4 @@ def convert_logger(loglev=None):
                         loglev, loglev_default)
         loglev = loglev_default
 
-    return loglev
\ No newline at end of file
+    return loglev
diff --git a/smmregrid/regrid.py b/smmregrid/regrid.py
index f176f6f..3eb5fb0 100644
--- a/smmregrid/regrid.py
+++ b/smmregrid/regrid.py
@@ -14,14 +14,15 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-# Modified slightly by Jost von Hardenberg <jost.hardenberg@polito.it>
+# Refactored  by Jost von Hardenberg <jost.hardenberg@polito.it>
+# and Paolo Davini <paolo.davini@cnr.it>
 
 """Dask-aware regridding
 
 To apply a regridding you will need a set of weights mapping from the source
 grid to the target grid.
 
-Regridding weights can be generated online using CDO (:func:`cdo_generate_weights`),
+Regridding weights can be generated online using CDO (:class:`CdoGenerate`),
 or offline by calling these programs externally (this is recommended especially
 for large grids).
 
@@ -34,319 +35,306 @@
 
 
 import math
-import sys
+import os
 import xarray
 import numpy
 import dask.array
 from .dimension import remove_degenerate_axes
-from .cdo_weights import cdo_generate_weights
-from .util import find_vert_coords
+from .cdogenerate import CdoGenerate
 from .weights import compute_weights_matrix3d, compute_weights_matrix, mask_weights, check_mask
 from .log import setup_logger
-
-# default spatial dimensions and vertical coordinates
-default_space_dims = ['i', 'j', 'x', 'y', 'lon', 'lat', 'longitude', 'latitude',
-                      'cell', 'cells', 'ncells', 'values', 'value', 'nod2', 'pix', 'elem']
-
-
-def apply_weights(source_data, weights, weights_matrix=None,
-                  masked=True, space_dims=None, loglevel='WARNING'):
-    """
-    Apply the CDO weights ``weights`` to ``source_data``, performing a regridding operation
-
-    Args:
-        source_data (xarray.DataArray): Source dataset
-        weights (xarray.DataArray): CDO weights information
-        masked (bool): if the DataArray is masked
-        space_dims (list): dimensions on which the interpolation has to be done (e.g.['lon', 'lat'])
-
-    Returns:
-        xarray.DataArray: Regridded version of the source dataset
-    """
-
-    loggy = setup_logger(level=loglevel, name='smmregrid.apply_weights')
-
-    # Understand immediately if we need to return something or not
-    # This is done if we have bounds variables
-    if any(substring in source_data.name for substring in ["bnds", "bounds", "vertices"]):
-
-        # we keep time bounds, and we ignore all the rest
-        if 'time' in source_data.name:
-            loggy.info('%s will not be interpolated in the output', source_data.name)
-            return source_data
-        else:
-            loggy.info('%s will be excluded from the output', source_data.name)
-            return xarray.DataArray(data=None)
-
-    # Alias the weights dataset from CDO
-    w = weights
-
-    # The weights file contains a sparse matrix, that we need to multiply the
-    # source data's horizontal grid with to get the regridded data.
-    #
-    # A bit of messing about with `.stack()` is needed in order to get the
-    # dimensions to conform - the horizontal grid needs to be converted to a 1d
-    # array, multiplied by the weights matrix, then unstacked back into a 2d
-    # array
-
-    # CDO style weights
-    # src_address = w.src_address - 1
-    # dst_address = w.dst_address - 1
-    # remap_matrix = w.remap_matrix[:, 0]
-    # w_shape = (w.sizes["src_grid_size"], w.sizes["dst_grid_size"])
-
-    dst_grid_shape = w.dst_grid_dims.values
-    dst_grid_center_lat = w.dst_grid_center_lat.data.reshape(
-        # dst_grid_shape[::-1], order="C"
-        dst_grid_shape[::-1]
-    )
-    dst_grid_center_lon = w.dst_grid_center_lon.data.reshape(
-        # dst_grid_shape[::-1], order="C"
-        dst_grid_shape[::-1]
-    )
-
-    dst_mask = w.dst_grid_imask
-    # src_mask = w.src_grid_imask
-
-    axis_scale = 180.0 / math.pi  # Weight lat/lon in radians
-
-    # Dimension on which we can produce the interpolation
-    if space_dims is None:
-        space_dims = default_space_dims
-
-    if not any(x in source_data.dims for x in space_dims):
-        loggy.error("None of dimensions on which we can interpolate is found in the DataArray. Does your DataArray include any of these?")
-        loggy.error(space_dims)
-        loggy.error('smmregrid can identify only %s', source_data.dims)
-        raise KeyError('Dimensions mismatch')
-
-    # Find dimensions to keep
-    nd = sum([(d not in space_dims) for d in source_data.dims])
-
-    kept_shape = list(source_data.shape[0:nd])
-    kept_dims = list(source_data.dims[0:nd])
-    loggy.info('Dimension kept: %s', kept_dims)
-
-    if weights_matrix is None:
-        weights_matrix = compute_weights_matrix(weights)
-
-    # Remove the spatial axes, apply the weights, add the spatial axes back
-    source_array = source_data.data
-    if isinstance(source_array, dask.array.Array):
-        source_array = dask.array.reshape(source_array, kept_shape + [-1])
-    else:
-        source_array = numpy.reshape(source_array, kept_shape + [-1])
-
-    # Handle input mask
-    dask.array.ma.set_fill_value(source_array, 1e20)
-    source_array = dask.array.ma.fix_invalid(source_array)
-    source_array = dask.array.ma.filled(source_array)
-
-    target_dask = dask.array.tensordot(source_array, weights_matrix, axes=1)
-
-    # define and compute the new mask
-    if masked:
-
-        # target mask is loaded from above
-        target_mask = dst_mask.data
-
-        # broadcast the mask on all the remaining dimensions
-        target_mask = numpy.broadcast_to(
-            target_mask.reshape([1 for d in kept_shape] + [-1]), target_dask.shape
-        )
-
-        # apply the mask
-        target_dask = dask.array.where(target_mask != 0.0, target_dask, numpy.nan)
-
-    # after the tensordot, bring the NaN back in
-    # Use greater than 1e19 to avoid numerical noise from interpolation.
-    target_dask = xarray.where(target_dask > 1e19, numpy.nan, target_dask)
-
-    # reshape the target DataArray
-    target_dask = dask.array.reshape(
-        target_dask, kept_shape + [dst_grid_shape[1], dst_grid_shape[0]]
-    )
-
-    # Create a new DataArray for the output
-    target_da = xarray.DataArray(
-        target_dask,
-        dims=kept_dims + ["i", "j"],
-        coords={
-            k: v
-            for k, v in source_data.coords.items()
-            if set(v.dims).issubset(kept_dims)
-        },
-        name=source_data.name,
-    )
-
-    target_da.coords["lat"] = xarray.DataArray(dst_grid_center_lat, dims=["i", "j"])
-    target_da.coords["lon"] = xarray.DataArray(dst_grid_center_lon, dims=["i", "j"])
-
-    # Clean up coordinates
-    target_da.coords["lat"] = remove_degenerate_axes(target_da.lat)
-    target_da.coords["lon"] = remove_degenerate_axes(target_da.lon)
-
-    # Convert to degrees if needed, rounding to avoid numerical errors
-    target_da.coords["lat"] = numpy.round(target_da.lat * axis_scale, 10)
-    target_da.coords["lon"] = numpy.round(target_da.lon * axis_scale, 10)
-
-    # If a regular grid drop the 'i' and 'j' dimensions
-    if target_da.coords["lat"].ndim == 1 and target_da.coords["lon"].ndim == 1:
-        target_da = target_da.swap_dims({"i": "lat", "j": "lon"})
-
-    # Add metadata to the coordinates
-    target_da.coords["lat"].attrs["units"] = "degrees_north"
-    target_da.coords["lat"].attrs["standard_name"] = "latitude"
-    target_da.coords["lon"].attrs["units"] = "degrees_east"
-    target_da.coords["lon"].attrs["standard_name"] = "longitude"
-    target_da.coords["lon"].attrs["axis"] = "X"
-    target_da.coords["lat"].attrs["axis"] = "Y"
-
-    # Copy attributes from the original
-    target_da.attrs = source_data.attrs
-
-    # Clean CDI gridtype (which can lead to issues with CDO interpretation)
-    if target_da.attrs.get('CDI_grid_type'):
-        del target_da.attrs['CDI_grid_type']
-
-    # Now rename to the original coordinate names
-    # target_da = target_da.rename({"lat": source_lat.name, "lon": source_lon.name})
-
-    return target_da
+from .gridinspector import GridInspector
+from .util import deprecated_argument
 
 
 class Regridder(object):
-    """Set up the regridding operation
+    """Main smmregrid regridding class"""
 
-    Supply either both ``source_grid`` and ``dest_grid`` or just ``weights``.
-
-    For large grids you may wish to pre-calculate the weights using
-    CDO, if not supplied ``weights`` will be calculated from
-    ``source_grid`` and ``dest_grid`` using CDO's genbil function.
+    def __init__(self, source_grid=None, target_grid=None, weights=None,
+                 method='con', transpose=True, vert_coord=None, vertical_dim=None,
+                 space_dims=None, horizontal_dims=None, cdo_extra=None, cdo_options=None,
+                 cdo='cdo', loglevel='WARNING'):
+        """
+        Initialize the Regridder for performing regridding operations.
 
-    Weights may be pre-computed by an external program, or created using
-    :func:`cdo_generate_weights`.
+        This class allows for regridding between two grids, either by providing
+        source and target grids or by using pre-computed weights. For large grids,
+        it is recommended to pre-calculate weights using CDO. If weights are not
+        supplied, they will be calculated using CDO's `genbil` function based on
+        the provided source and target grids.
 
-    Args:
-        source_grid (:class:`xarray.DataArray`): Source grid / sample dataset
-        target_grid (:class:`xarray.DataArray`): Target grid / sample dataset
-        weights (:class:`xarray.Dataset`): Pre-computed interpolation weights
-        vert_coord (str): Name of the vertical coordinate.
-                          If provided, 3D weights are generated (default: None)
-        level_idx (str): Prefix of helper vertical coordinate with original level indices.
-                         If provided, 3D weights are selected from those levels (default: "idx_")
-        method (str): Method to use for interpolation (default: 'con')
-        space_dims (list): list of dimensions to interpolate (default: None)
-        transpose (bool): transpose the output so that the vertical coordinate is
-                          just before other spatial coords (dafault: True)
-    """
+        Pre-computed weights can be generated externally or by using
+        :class:`CdoGenerate()`.
 
-    def __init__(self, source_grid=None, target_grid=None, weights=None,
-                 method='con', space_dims=None, vert_coord=None, transpose=True,
-                 cdo='cdo', level_idx="idx_", loglevel='WARNING'):
+        Args:
+            source_grid (xarray.DataArray or str): Source grid dataset or file path.
+            target_grid (xarray.DataArray): Target grid dataset for regridding.
+            weights (xarray.Dataset): Pre-computed interpolation weights.
+            vertical_dim (str): Name of the vertical coordinate for 
+                                3D weights generation (default: None).
+            horizontal_dims (list): List of spatial dimensions to 
+                                    interpolate (default: None).
+            method (str): Interpolation method to use (default: 'con').
+            transpose (bool): If True, transpose the output such that 
+                              the vertical coordinate is placed just before the 
+                              other spatial coordinates (default: True).
+            cdo (str): Path to the CDO executable (default: 'cdo').
+            cdo_extra (str): Extra command to be passed to CDO
+            cdo_options(str): Extra options to be passed to CDO
+            loglevel (str): Logging level for the operation (default: 'WARNING').
+
+        Raises:
+            ValueError: If neither weights nor source and target grids are provided.
+            FileNotFoundError: If the specified source grid file does not exist.
+            KeyError: If no grid types are found in the data.
+
+        Warnings:
+            DeprecationWarning: If deprecated arguments 'vert_coord' or 'space_dims' are used.
+        """
 
         if (source_grid is None or target_grid is None) and (weights is None):
             raise ValueError(
                 "Either weights or source_grid/target_grid must be supplied"
             )
-        
-        self.loggy = setup_logger(level=loglevel, name='smmregrid.regrid')
-        self.loglevel = loglevel
 
+        vertical_dim = deprecated_argument(vert_coord, vertical_dim, 'vert_coord', 'vertical_dim')
+        horizontal_dims = deprecated_argument(space_dims, horizontal_dims, 'space_dims', 'horizontal_dims')
+       
+        # set up logger
+        self.loggy = setup_logger(level=loglevel, name='smmregrid.Regrid')
+        self.loglevel = loglevel
         self.transpose = transpose
+        vertical_dim = [vertical_dim] if isinstance(vertical_dim, str) else vertical_dim
+        horizontal_dims = [horizontal_dims] if isinstance(horizontal_dims, str) else horizontal_dims
+        if vertical_dim:
+            self.loggy.info('Forcing vertical_dim from input: expecting a single-gridtype dataset')
+        if horizontal_dims:
+            self.loggy.info('Forcing horizontal_dim from input: expecting a single-gridtype dataset')
+        self.extra_dims = {'vertical': vertical_dim, 'horizontal': horizontal_dims}
 
         # Is there already a weights file?
         if weights is not None:
-            if not isinstance(weights, xarray.Dataset):
-                self.weights = xarray.open_mfdataset(weights)
+            self.loggy.info('Init from weights selected!')
+            self.init_mode = 'weights'
+            self.grids = self._gridtype_from_weights(weights)
+        else:
+            self.loggy.info('Init from grids selected!')
+            self.init_mode = 'grids'
+            if isinstance(source_grid, str):
+                if os.path.isfile(source_grid):
+                    source_grid_array = xarray.open_dataset(source_grid)
+                else:
+                    raise FileNotFoundError(f'Cannot find grid file {source_grid}')
             else:
-                self.weights = weights
+                source_grid_array = source_grid
+
+            self.grids = self._gridtype_from_data(source_grid_array)
 
-            if not vert_coord:
-                self.vert_coord = find_vert_coords(self.weights)
+            len_grids = len(self.grids)
+            if len_grids == 0:
+                raise KeyError('Cannot find any gridtype in your data, aborting!')
+            if len_grids == 1:
+                self.loggy.info('One gridtype found! Standard procedure')
             else:
-                self.vert_coord = vert_coord
-        else:
+                self.loggy.warning('%s gridtypes found! We are in uncharted territory!', len_grids)
+
+            for index, gridtype  in enumerate(self.grids):
+                self.loggy.info('Processing grid number %s', index)
+                self.loggy.debug('Processing grids %s', gridtype.dims)
+                self.loggy.debug('Horizontal dimensions are %s', gridtype.horizontal_dims)
+                self.loggy.debug('Vertical dimension is %s', gridtype.vertical_dim)
+                self.loggy.debug('Other dimensions are %s', gridtype.other_dims)
+
+                # always prefer to pass filename (i.e. source_grid) when possible to CdoGenerate()
+                # this will limit errors from xarray and speed up CDO itself
+                # it wil work only for single-gridtype dataset
+                if isinstance(source_grid, str) and len_grids == 1:
+                    source_grid_array_to_cdo = source_grid
+                else:
+                    # when feeding from xarray, select the variable and its bounds
+                    if isinstance(source_grid_array, xarray.Dataset):
+                        stored_vars = [list(gridtype.variables.keys())[0]] + gridtype.bounds
+                        self.loggy.debug('Storing variables %s', stored_vars)
+                        source_grid_array_to_cdo = source_grid_array[stored_vars]
+                    else:
+                        source_grid_array_to_cdo = source_grid_array
+                    
+                    if gridtype.time_dims:
+                        self.loggy.debug('Selecting only first time step for dimension %s', gridtype.time_dims[0])
+                        source_grid_array_to_cdo = source_grid_array_to_cdo.isel({gridtype.time_dims[0]: 0})
+
+                generator = CdoGenerate(source_grid_array_to_cdo, target_grid,
+                                               cdo=cdo, cdo_options=cdo_options,
+                                               cdo_extra=cdo_extra, loglevel=loglevel)
+                gridtype.weights = generator.weights(method=method,
+                                                     vertical_dim=gridtype.vertical_dim)
+
+        for gridtype in self.grids:
+            if gridtype.vertical_dim:
+                gridtype.weights_matrix = compute_weights_matrix3d(gridtype.weights,
+                                                                   gridtype.vertical_dim)
+            else:
+                gridtype.weights_matrix = compute_weights_matrix(gridtype.weights)
 
-            # Check if there is a vertical coordinate for 3d oceanic data
-            if not vert_coord:
-                self.vert_coord = find_vert_coords(source_grid)
+            # this section is used to create a target mask initializing the CDO weights (both 2d and 3d)
+            # has a destination mask been precomputed?
+            if "dst_grid_masked" in gridtype.weights.variables:
+                gridtype.masked = gridtype.weights.dst_grid_masked.data
             else:
-                self.vert_coord = vert_coord
+                # compute the destination mask now
+                gridtype.weights = mask_weights(gridtype.weights, gridtype.weights_matrix, gridtype.vertical_dim)
+                gridtype.masked = check_mask(gridtype.weights, gridtype.vertical_dim)
+
+    def _gridtype_from_weights(self, weights):
+        """
+        Initialize the gridtype reading from weights
+        """
 
-            # Generate the weights with CDO
-            self.weights = cdo_generate_weights(source_grid, target_grid, method=method,
-                                                vert_coord=self.vert_coord, cdo=cdo)
+        self.loggy.info('Precomputed weights support so far single-gridtype datasets')
 
-        if self.vert_coord:
-            self.weights_matrix = compute_weights_matrix3d(self.weights, self.vert_coord)
-        else:
-            self.weights_matrix = compute_weights_matrix(self.weights)
+        if not isinstance(weights, xarray.Dataset):
+            weights = xarray.open_mfdataset(weights)
 
-        # this section is used to create a target mask initializing the CDO weights (both 2d and 3d)
-        if "dst_grid_masked" in self.weights.variables:  # has a destination mask been precomputed?
-            self.masked = self.weights.dst_grid_masked.data  # ok, let's use it
-        else:
-            # compute the destination mask now
-            self.weights = mask_weights(self.weights, self.weights_matrix, self.vert_coord)
-            self.masked = check_mask(self.weights, self.vert_coord)
+        grid_info = GridInspector(weights, cdo_weights=True, extra_dims=self.extra_dims,
+                                  clean=False, loglevel=self.loglevel)
+        gridtype = grid_info.get_grid_info()
+
+        # the vertical dimension has to show up into the extra dimensions
+        # to cover the case that it is not a standard dimension: possibly better implementation available
+        self.extra_dims['vertical'] = [gridtype[0].vertical_dim]
 
-        self.space_dims = space_dims
-        self.level_idx = level_idx
+        return gridtype
+
+    def _gridtype_from_data(self, source_grid_array):
+        """
+        Initialize the gridtype reading from source_data
+        """
+        grid_info = GridInspector(source_grid_array, extra_dims=self.extra_dims,
+                                  clean=True, loglevel=self.loglevel)
+        return grid_info.get_grid_info()
 
     def regrid(self, source_data):
-        """Regrid ``source_data`` to match the target grid
+        """
+        Regrid the provided source data to match the target grid.
+
+        This method applies the regridding process to either a single data array or
+        a dataset containing multiple data variables.
 
         Args:
-            source_data (:class:`xarray.DataArray` or xarray.Dataset): Source
-            variable
+            source_data (xarray.DataArray or xarray.Dataset): The source data to be regridded.
 
         Returns:
-            :class:`xarray.DataArray` or xarray.Dataset with a regridded
-            version of the source variable
+            xarray.DataArray or xarray.Dataset: The regridded data, matching the target grid.
+
+        Raises:
+            TypeError: If the provided source data is neither an xarray DataArray nor Dataset.
         """
 
+        # apply the regridder on each DataArray
         if isinstance(source_data, xarray.Dataset):
 
-            # apply the regridder on each DataArray
-            out = source_data.map(self.regrid_array, keep_attrs=True)
+            # extra call to GridInspector for the raise with multiple grids
+            grid_inspect = GridInspector(source_data,
+                                     extra_dims=self.extra_dims,
+                                     loglevel=self.loglevel)
+            datagrids = grid_inspect.get_grid_info()
+            if len(datagrids)>1 and self.init_mode == 'weights':
+                raise ValueError(f'Cannot process data with {len(datagrids)} GridType initializing from weights')
+
+            # map on multiple dataarray
+            out = source_data.map(self.regrid_array, keep_attrs=False)
 
             # clean from degenerated variables
             degen_vars = [var for var in out.data_vars if out[var].dims == ()]
             return out.drop_vars(degen_vars)
 
-        elif isinstance(source_data, xarray.DataArray):
-
+        if isinstance(source_data, xarray.DataArray):
             return self.regrid_array(source_data)
 
-        else:
-            sys.exit('The object provided is not a Xarray object!')
+        raise TypeError('The object provided is not a Xarray object!')
 
     def regrid_array(self, source_data):
-        """Regridding selection through 2d and 3d arrays"""
+        """
+        Perform regridding on a single 2D or 3D array.
 
-        if (self.vert_coord and self.vert_coord in source_data.coords):
-            # if this is a 3D we specified the vertical coord and it has it
-            return self.regrid3d(source_data)
-        else:
-            return self.regrid2d(source_data)
+        This method decides whether to call the 2D or 3D regridding function based
+        on the dimensions of the source data.
+
+        Args:
+            source_data (xarray.DataArray): The source data array to regrid.
+
+        Returns:
+            xarray.DataArray: The regridded array.
+
+        Raises:
+            ValueError: If the input data does not match expected dimensions.
+        """
+        
+        self.loggy.debug('Getting GridType from source_data')
+        grid_inspect = GridInspector(source_data,
+                                     extra_dims=self.extra_dims,
+                                     loglevel=self.loglevel)
+        datagrids = grid_inspect.get_grid_info()
+
+        for datagridtype in datagrids:
+            if datagridtype.vertical_dim:
+                # if this is a 3D we specified the vertical coord and it has it
+                return self.regrid3d(source_data, datagridtype)
+            # 2d case
+            return self.regrid2d(source_data, datagridtype)
+
+    def _get_gridtype(self, datagridtype):
+
+        # special case for CDO weights without any dimensional information
+        # we derived this from the regridded data and we use it as it is
+        if self.init_mode == 'weights':
+            self.loggy.info('Assuming gridtype from data to be the same from weights')
+            self.grids[0].dims = datagridtype.dims
+            self.grids[0].horizontal_dims = datagridtype.horizontal_dims
+            self.grids[0].other_dims = datagridtype.other_dims
+
+        # match the grid
+        gridtype = next((grid for grid in self.grids if grid == datagridtype), None)
+
+        return gridtype
+
+    def regrid3d(self, source_data, datagridtype):
+        """
+        Regrid a 3D source data array to match the target grid.
 
-    def regrid3d(self, source_data):
-        """Regrid ``source_data`` to match the target grid - 3D version
+        This method applies the necessary weights and handles vertical coordinates.
 
         Args:
-            source_data (:class:`xarray.DataArray`): Source
-            variable
+            source_data (xarray.DataArray): The 3D source data to be regridded.
+            datagridtype: The grid type information for the source data.
 
         Returns:
-            :class:`xarray.DataArray` with a regridded
-            version of the source variable
+            xarray.DataArray: The regridded 3D data.
+
+        Raises:
+            ValueError: If there are issues with transposing the output dimensions.
         """
 
         self.loggy.debug('3D DataArray access: variable is %s', source_data.name)
 
+        # get the gridtype from class that matches the data
+        gridtype = self._get_gridtype(datagridtype)
+
+        # if the current grid is not available in the Regrid gridtype
+        if gridtype is None:
+            self.loggy.info('%s will be excluded from the output', source_data.name)
+            return xarray.DataArray(data=None)
+
+        # select the gridtype to be used
+        vertical_dim = gridtype.vertical_dim
+        weights = gridtype.weights
+        weights_matrix = gridtype.weights_matrix
+        masked = gridtype.masked
+        level_index = gridtype.level_index
+        horizontal_dims = gridtype.horizontal_dims
+
         # CDO 2.2.0 fix
-        if "numLinks" in self.weights.dims:
+        if "numLinks" in weights.dims:
             links_dim = "numLinks"
         else:
             links_dim = "num_links"
@@ -357,62 +345,258 @@ def regrid3d(self, source_data):
             return source_data
 
         # If a special additional coordinate is present pick correct levels from weights
-        coord = next((coord for coord in source_data.coords if coord.startswith(self.level_idx)), None)
-        if coord:  # if a coordinate starting with level_idx is found
+        coord = next((coord for coord in source_data.coords if coord.startswith(level_index)), None)
+        if coord:  # if a coordinate starting with level_index is found
             levlist = source_data.coords[coord].values.tolist()
             levlist = [levlist] if numpy.isscalar(levlist) else levlist
         else:
-            levlist = list(range(0, source_data.coords[self.vert_coord].values.size))
+            levlist = list(range(0, source_data.coords[vertical_dim].values.size))
 
         data3d_list = []
         for lev, levidx in enumerate(levlist):
-            xa = source_data.isel(**{self.vert_coord: lev})
-            wa = self.weights.isel(**{self.vert_coord: levidx})
+            self.loggy.debug('Processing vertical level %s - level_index %s', lev, levidx)
+            xa = source_data.isel(**{vertical_dim: lev})
+            wa = weights.isel(**{vertical_dim: levidx})
             nl = wa.link_length.values
             wa = wa.isel(**{links_dim: slice(0, nl)})
-            wm = self.weights_matrix[levidx]
-            mm = self.masked[levidx]
-            data3d_list.append(apply_weights(
+            wm = weights_matrix[levidx]
+            mm = masked[levidx]
+            data3d_list.append(self.apply_weights(
                 xa, wa, weights_matrix=wm,
-                masked=mm, space_dims=self.space_dims,
-                loglevel=self.loglevel)
+                masked=mm, horizontal_dims=horizontal_dims)
             )
-        data3d = xarray.concat(data3d_list, dim=self.vert_coord)
+        data3d = xarray.concat(data3d_list, dim=vertical_dim)
+
+        # get dimensional info on target grid. TODO: can be moved at the init?
+        target_gridtypes = GridInspector(data3d, clean=True, loglevel=self.loglevel).get_grid_info()
+        target_horizontal_dims = target_gridtypes[0].horizontal_dims
 
         if self.transpose:
-            # Make sure that the vertical dimension is just before the spatial ones
-            if self.space_dims is None:
-                space_dims = default_space_dims
-            else:
-                space_dims = self.space_dims
             dims = list(data3d.dims)
-            index = min([i for i, s in enumerate(dims) if s in space_dims])
+            index = min([i for i, s in enumerate(dims) if s in target_horizontal_dims])
             dimst = dims[1:index] + [dims[0]] + dims[index:]
             data3d = data3d.transpose(*dimst)
 
             return data3d
         else:
-            raise ValueError('Cannot process this source data, are you sure it is an xarray?')
+            raise ValueError(f'Cannot transpose output dimensions {data3d.dims} over {target_horizontal_dims}')
+
+    def regrid2d(self, source_data, datagridtype):
+        """
+        Regrid a 2D source data array to match the target grid.
 
-    def regrid2d(self, source_data):
-        """Regrid ``source_data`` to match the target grid, 2D version
+        This method applies the necessary weights for 2D data regridding.
 
         Args:
-            source_data (:class:`xarray.DataArray`): Source
-            variable
+            source_data (xarray.DataArray): The 2D source data to be regridded.
+            datagridtype: The grid type information for the source data.
 
         Returns:
-            :class:`xarray.DataArray` with a regridded
-            version of the source variable
+            xarray.DataArray: The regridded 2D data.
         """
         self.loggy.debug('2D DataArray access: variables is %s', source_data.name)
-        return apply_weights(
-            source_data, self.weights, weights_matrix=self.weights_matrix,
-            masked=self.masked, space_dims=self.space_dims, loglevel=self.loglevel
+
+        # get the gridtype from class that matches the data
+        gridtype = self._get_gridtype(datagridtype)
+
+        if gridtype is None:
+            self.loggy.info('%s will be excluded from the output', source_data.name)
+            return xarray.DataArray(data=None)
+
+        return self.apply_weights(
+            source_data,
+            gridtype.weights,
+            weights_matrix=gridtype.weights_matrix,
+            masked=gridtype.masked,
+            horizontal_dims=gridtype.horizontal_dims)
+
+    def apply_weights(self, source_data, weights, weights_matrix=None,
+                      masked=True, horizontal_dims=None):
+        """
+        Apply CDO weights to the source data, performing the regridding operation.
+
+        This method multiplies the source data with the weights matrix to produce
+        the regridded output.
+
+        Args:
+            source_data (xarray.DataArray): The source data to be regridded.
+            weights (xarray.DataArray): The CDO weights for interpolation.
+            weights_matrix (xarray.DataArray, optional): Pre-computed weights matrix (default: None).
+            masked (bool): Indicates if the DataArray is masked (default: True).
+            horizontal_dims (list): List of dimensions for interpolation (e.g., ['lon', 'lat']).
+
+        Returns:
+            xarray.DataArray: The regridded version of the source dataset.
+
+        Raises:
+            KeyError: If none of the specified horizontal dimensions are found in the DataArray.
+        """
+
+        # Understand immediately if we need to return something or not
+        # This is done if we have bounds variables
+        if any(substring in source_data.name for substring in ["bnds", "bounds", "vertices"]):
+
+            # we keep time bounds, and we ignore all the rest
+            if 'time' in source_data.name:
+                self.loggy.info('%s will not be interpolated in the output', source_data.name)
+                return source_data
+            else:
+                self.loggy.info('%s will be excluded from the output', source_data.name)
+                return xarray.DataArray(data=None)
+
+        # Alias the weights dataset from CDO
+        w = weights
+
+        # The weights file contains a sparse matrix, that we need to multiply the
+        # source data's horizontal grid with to get the regridded data.
+        #
+        # A bit of messing about with `.stack()` is needed in order to get the
+        # dimensions to conform - the horizontal grid needs to be converted to a 1d
+        # array, multiplied by the weights matrix, then unstacked back into a 2d
+        # array
+
+        # CDO style weights
+        # src_address = w.src_address - 1
+        # dst_address = w.dst_address - 1
+        # remap_matrix = w.remap_matrix[:, 0]
+        # w_shape = (w.sizes["src_grid_size"], w.sizes["dst_grid_size"])
+
+        # src_grid_rank = w.src_grid_rank
+        dst_grid_rank = w.dst_grid_rank
+
+        src_cdo_grid = w.attrs['source_grid']
+        dst_cdo_grid = w.attrs['dest_grid']
+        self.loggy.info('Interpolating from CDO %s to CDO %s', src_cdo_grid, dst_cdo_grid)
+
+        dst_grid_shape = w.dst_grid_dims.values
+        dst_grid_center_lat = w.dst_grid_center_lat.data.reshape(
+            dst_grid_shape[::-1]
         )
+        dst_grid_center_lon = w.dst_grid_center_lon.data.reshape(
+            dst_grid_shape[::-1]
+        )
+
+        dst_mask = w.dst_grid_imask
+        # src_mask = w.src_grid_imask
+
+        axis_scale = 180.0 / math.pi  # Weight lat/lon in radians
+
+        if not any(x in source_data.dims for x in horizontal_dims):
+            self.loggy.error(
+                "None of dimensions on which we can interpolate is found in the DataArray. Does your DataArray include any of these?")
+            self.loggy.error(horizontal_dims)
+            self.loggy.error('smmregrid can identify only %s', source_data.dims)
+            raise KeyError('Dimensions mismatch')
+
+        self.loggy.info('Regridding from %s to %s', source_data.shape, dst_grid_shape)
+
+        # Find dimensions to keep
+        nd = sum([(d not in horizontal_dims) for d in source_data.dims])
+
+        kept_shape = list(source_data.shape[0:nd])
+        kept_dims = list(source_data.dims[0:nd])
+        self.loggy.debug('Dimension to be ignored: %s', kept_dims)
+        if weights_matrix is None:
+            weights_matrix = compute_weights_matrix(weights)
+
+        # Remove the spatial axes, apply the weights, add the spatial axes back
+        source_array = source_data.data
+        if isinstance(source_array, dask.array.Array):
+            source_array = dask.array.reshape(source_array, kept_shape + [-1])
+        else:
+            source_array = numpy.reshape(source_array, kept_shape + [-1])
+        self.loggy.debug('Source array after reshape is: %s', source_array.shape)
+
+        # Handle input mask
+        dask.array.ma.set_fill_value(source_array, 1e20)
+        source_array = dask.array.ma.fix_invalid(source_array)
+        source_array = dask.array.ma.filled(source_array)
+
+        self.loggy.debug('Tensordot!')
+        target_dask = dask.array.tensordot(source_array, weights_matrix, axes=1)
+
+        # define and compute the new mask
+        if masked:
+
+            # target mask is loaded from above
+            target_mask = dst_mask.data
+
+            # broadcast the mask on all the remaining dimensions
+            target_mask = numpy.broadcast_to(
+                target_mask.reshape([1 for d in kept_shape] + [-1]), target_dask.shape
+            )
+            self.loggy.debug('Reshaped mask with %s', target_mask.shape)
+
+            # apply the mask
+            target_dask = dask.array.where(target_mask != 0.0, target_dask, numpy.nan)
+
+        # after the tensordot, bring the NaN back in
+        # Use greater than 1e19 to avoid numerical noise from interpolation.
+        target_dask = xarray.where(target_dask > 1e19, numpy.nan, target_dask)
+
+        if len(dst_grid_rank) == 2:
+            tgt_shape = [dst_grid_shape[1], dst_grid_shape[0]]
+            tgt_dims = ["i", "j"]
+        elif len(dst_grid_rank) == 1:
+            tgt_shape = [dst_grid_shape[0]]
+            tgt_dims = ['cell']
+        else:
+            raise ValueError('Unknown dimensional target grid')
+
+        # reshape the target DataArray
+        target_dask = dask.array.reshape(
+            target_dask, kept_shape + tgt_shape
+        )
+
+        # Create a new DataArray for the output
+        target_da = xarray.DataArray(
+            target_dask,
+            dims=kept_dims + tgt_dims,
+            coords={
+                k: v
+                for k, v in source_data.coords.items()
+                if set(v.dims).issubset(kept_dims)
+            },
+            name=source_data.name,
+        )
+
+        target_da.coords["lat"] = xarray.DataArray(dst_grid_center_lat, dims=tgt_dims)
+        target_da.coords["lon"] = xarray.DataArray(dst_grid_center_lon, dims=tgt_dims)
+
+        # Clean up coordinates
+        target_da.coords["lat"] = remove_degenerate_axes(target_da.lat)
+        target_da.coords["lon"] = remove_degenerate_axes(target_da.lon)
+
+        # Convert to degrees if needed, rounding to avoid numerical errors
+        target_da.coords["lat"] = numpy.round(target_da.lat * axis_scale, 10)
+        target_da.coords["lon"] = numpy.round(target_da.lon * axis_scale, 10)
+
+        # If a regular grid drop the 'i' and 'j' dimensions
+        if tgt_dims == ['i', 'j'] and target_da.coords["lat"].ndim == 1 and target_da.coords["lon"].ndim == 1:
+            target_da = target_da.swap_dims({"i": "lat", "j": "lon"})
+
+        # Add metadata to the coordinates
+        target_da.coords["lat"].attrs["units"] = "degrees_north"
+        target_da.coords["lat"].attrs["standard_name"] = "latitude"
+        target_da.coords["lon"].attrs["units"] = "degrees_east"
+        target_da.coords["lon"].attrs["standard_name"] = "longitude"
+        target_da.coords["lon"].attrs["axis"] = "X"
+        target_da.coords["lat"].attrs["axis"] = "Y"
+
+        # Copy attributes from the original
+        target_da.attrs = source_data.attrs
+
+        # Clean CDI gridtype (which can lead to issues with CDO interpretation)
+        if target_da.attrs.get('CDI_grid_type'):
+            del target_da.attrs['CDI_grid_type']
+
+        # Now rename to the original coordinate names
+        # target_da = target_da.rename({"lat": source_lat.name, "lon": source_lon.name})
+
+        return target_da
 
 
-def regrid(source_data, target_grid=None, weights=None, vert_coord=None, transpose=True, cdo='cdo'):
+def regrid(source_data, target_grid=None, weights=None, transpose=True, cdo='cdo'):
     """
     A simple regrid. Inefficient if you are regridding more than one dataset
     to the target grid because it re-generates the weights each time you call
@@ -423,8 +607,6 @@ def regrid(source_data, target_grid=None, weights=None, vert_coord=None, transpo
     Args:
         source_data (:class:`xarray.DataArray`): Source variable
         target_grid (:class:`coecms.grid.Grid` or :class:`xarray.DataArray`): Target grid / sample variable
-        vert_coord (str): Name of the vertical coordinate.
-                          If provided, 3D weights are generated (default: None)
         weights (:class:`xarray.Dataset`): Pre-computed interpolation weights
         transpose (bool): If True, transpose the output so that the vertical
                           coordinate is just before the other spatial coordinates (default: True)
@@ -434,8 +616,7 @@ def regrid(source_data, target_grid=None, weights=None, vert_coord=None, transpo
         :class:`xarray.DataArray` with a regridded version of the source variable
     """
 
-    regridder = Regridder(source_data, target_grid=target_grid, weights=weights,
-                          vert_coord=vert_coord, cdo=cdo, transpose=transpose)
+    regridder = Regridder(source_data, target_grid=target_grid, weights=weights, cdo=cdo, transpose=transpose)
     return regridder.regrid(source_data)
 
 
diff --git a/smmregrid/util.py b/smmregrid/util.py
index c7fe534..fa6d64c 100644
--- a/smmregrid/util.py
+++ b/smmregrid/util.py
@@ -1,26 +1,64 @@
-"""Utils module"""
+"""Utility module with various support function"""
 
+import re
+import os
+import warnings
 import xarray
-import logging
 
-default_vert_coords = ['lev', 'nz1', 'nz', 'depth', 'depth_full', 'depth_half']
+# Define CDO regex patterns for each grid type (updated at CDO 2.4.4)
+# Define regex patterns for each grid type
+CDO_GRID_PATTERNS = {
+    "global_regular": re.compile(r"^global_\d+(\.\d+)?$"),                    
+    "regional_regular": re.compile(r"^dcw:[A-Z]{2,4}(?:_\d+(\.\d+)?)?$"),     
+    "zonal_latitudes": re.compile(r"^zonal_\d+(\.\d+)?$"),                    
+    "global_regular_NxM": re.compile(r"^r\d+x\d+$"),                          
+    "one_grid_point": re.compile(r"^lon=(-?\d+(\.\d+)?)/lat=(-?\d+(\.\d+)?)$"),
+    "gaussian_grid_F": re.compile(r"^F\d+$"),      
+    "gaussian_grid_n": re.compile(r"^n\d+$"),                            
+    "icosahedral_gme": re.compile(r"^gme\d+$"),                           
+    "healpix_grid": re.compile(r"^hp\d+(?:_(nested|ring))?$"),             
+    "healpix_zoom": re.compile(r"^hpz\d+$")                       
+}
 
-# set up logger
-loggy = logging.getLogger(__name__)
+def is_cdo_grid(grid_str):
+    """Check if the input string matches any CDO grid type."""
 
+    # Check if the string matches any of the grid patterns
+    for grid_type, pattern in CDO_GRID_PATTERNS.items():
+        if pattern.match(grid_str):
+            return True
 
-def find_vert_coords(xfield):
-    """
-    Find a vertical coordinate among defaults
-    Used to define if we need the 3d interpolation with adaptive mask
-    """
+    # Return False if no patterns match
+    return False
 
-    if isinstance(xfield, str):
-        xfield = xarray.open_dataset(xfield)
+def check_gridfile(filename):
+    """Check if a grid is a file, a cdo string or xarray object"""
 
-    vcoords = list(set(xfield.coords.keys()).intersection(default_vert_coords))
-    if len(vcoords) > 1:
-        raise ValueError('Multiple vertical coordinates in the same file are not yet supported')
-    if len(vcoords) == 0:
+    if filename is None:
         return None
-    return vcoords[0]
+    if isinstance(filename, (xarray.Dataset, xarray.DataArray)):
+        return "xarray"
+    if isinstance(filename, str):
+        if is_cdo_grid(filename):
+            return "grid"
+        if os.path.exists(filename):
+            return "file"
+        raise FileNotFoundError(f"Cannot find {filename} on disk")
+
+    raise TypeError(f'Unsuported format for {filename}')
+
+
+def deprecated_argument(old, new, oldname='var1', newname='var2'):
+    """Utility to provide warning in case of deprecated argument"""
+
+    # Check for deprecated 'old' argument
+    if old is not None:
+        warnings.warn(
+            f"{oldname} is deprecated and will be removed in future versions. "
+            f"Please use {newname} instead.",
+            DeprecationWarning
+        )
+        # If new is not provided, use the value from old
+        if new is None:
+            new = old
+    return new
diff --git a/smmregrid/weights.py b/smmregrid/weights.py
index 9637e9b..9e37b13 100644
--- a/smmregrid/weights.py
+++ b/smmregrid/weights.py
@@ -2,12 +2,9 @@
 
 import dask.array
 import sparse
-import logging
 
-loggy = logging.getLogger(__name__)
 
-
-def compute_weights_matrix3d(weights, vert_coord='lev'):
+def compute_weights_matrix3d(weights, vertical_dim='lev'):
     """
     Convert the weights from CDO to a list of numpy arrays
     """
@@ -19,9 +16,9 @@ def compute_weights_matrix3d(weights, vert_coord='lev'):
         links_dim = "num_links"
 
     sparse_weights = []
-    nvert = weights[vert_coord].values.size
+    nvert = weights[vertical_dim].values.size
     for i in range(0, nvert):
-        w = weights.loc[{vert_coord: i}]
+        w = weights.loc[{vertical_dim: i}]
         nl = w.link_length.values
         w = w.isel(**{links_dim: slice(0, nl)})
         sparse_weights.append(compute_weights_matrix(w))
@@ -60,16 +57,16 @@ def mask_tensordot(src_mask, weights_matrix):
     return target_mask
 
 
-def mask_weights(weights, weights_matrix, vert_coord=None):
+def mask_weights(weights, weights_matrix, vertical_dim=None):
     """This functions precompute the mask for the target interpolation
     Takes as input the weights from CDO and the precomputed weights matrix
     Return the target mask: handle the 3d case"""
 
     src_mask = weights.src_grid_imask
-    if vert_coord is not None:
-        for nlev in range(len(weights[vert_coord])):
-            mask = src_mask.loc[{vert_coord: nlev}].data
-            weights['dst_grid_imask'].loc[{vert_coord: nlev}] = mask_tensordot(mask, weights_matrix[nlev])
+    if vertical_dim is not None:
+        for nlev in range(len(weights[vertical_dim])):
+            mask = src_mask.loc[{vertical_dim: nlev}].data
+            weights['dst_grid_imask'].loc[{vertical_dim: nlev}] = mask_tensordot(mask, weights_matrix[nlev])
     else:
         mask = src_mask.data
         weights['dst_grid_imask'].data = mask_tensordot(mask, weights_matrix)
@@ -77,41 +74,17 @@ def mask_weights(weights, weights_matrix, vert_coord=None):
     return weights
 
 
-def check_mask(weights, vert_coord=None):
+def check_mask(weights, vertical_dim=None):
     """Check if the target mask is empty or full and
     return a bool to be passed to the regridder.
     Handle the 3d case (5x time faster)"""
 
     wdst = weights['dst_grid_imask']
-    if vert_coord is not None:
-        check = wdst.mean(dim=tuple(dim for dim in wdst.dims if dim != vert_coord))
+    if vertical_dim is not None:
+        check = wdst.mean(dim=tuple(dim for dim in wdst.dims if dim != vertical_dim))
         out = (check != 1).data.tolist()
     else:
         check = wdst.mean()
         out = (check != 1).data
 
     return out
-
-
-# def check_mask_old(weights, vert_coord=None):
-#     """Old version: check if the target mask is empty or full and
-#     return a bool to be passed to the regridder.
-#     Handle the 3d case"""
-
-#     if vert_coord is not None:
-#         result = []
-#         for nlev in range(len(weights[vert_coord])):
-#             w = weights['dst_grid_imask'].loc[{vert_coord: nlev}]
-#             v = w.sum() / len(w)
-#             if v == 1:
-#                 result.append(False)
-#             else:
-#                 result.append(True)
-#         return result
-#     else:
-#         w = weights['dst_grid_imask']
-#         v = w.sum() / len(w)
-#         if v == 1:
-#             return False
-#         else:
-#             return True
diff --git a/playground.ipynb b/speed-evaluation.ipynb
similarity index 52%
rename from playground.ipynb
rename to speed-evaluation.ipynb
index 1898b1c..4c14873 100644
--- a/playground.ipynb
+++ b/speed-evaluation.ipynb
@@ -5,23 +5,22 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Tests for SMM (without dask) versus CDO\n",
+    "# Speed tests for smmregrid (without dask) versus CDO\n",
     "\n",
-    "There are a few tests to check if the SMM approach is faster than the CDO one and if it is reliable in terms of output. Encouraging results below. Tested with both 2D and 3D vars, using DataArray and Datasets"
+    "There are a few tests to check if the SMM approach is faster than the CDO one and if it is reliable in terms of output. Tested with both 2D and 3D vars, using DataArray and Datasets.\n",
+    "More complex oceanic 3D data structures are not yet considered"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
     "import timeit\n",
     "import os\n",
-    "import numpy as np\n",
     "import xarray as xr\n",
     "from smmregrid import cdo_generate_weights, Regridder\n",
-    "from smmregrid.checker import check_cdo_regrid # this is a new function introduced to verify the output\n",
     "from cdo import Cdo\n",
     "import pandas as pd\n",
     "import copy\n",
@@ -31,8 +30,9 @@
     "\n",
     "# where and which the data are\n",
     "indir='tests/data'\n",
-    "filelist = ['tas-healpix2.nc', 'lsm-ifs.grb','onlytos-ipsl.nc','tas-ecearth.nc', \n",
+    "filelist = ['tas-healpix2.nc', 'onlytos-ipsl.nc','tas-ecearth.nc', \n",
     "            '2t-era5.nc','tos-fesom.nc', 'ua-ecearth.nc', 'mix-cesm.nc']#,'era5-mon.nc'] # the last is not available on github\n",
+    "#'lsm-ifs.grb'\n",
     "#filelist = ['tos-fesom.nc','onlytos-ipsl.nc','tas-ecearth.nc'] \n",
     "#filelist = ['tas-ecearth.nc']\n",
     "tfile = os.path.join(indir, 'r360x180.nc')\n",
@@ -62,109 +62,6 @@
     "    base['ua-ecearth.nc']['chunks'] = {'plev': 3}"
    ]
   },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Robustness test\n",
-    "\n",
-    "This is to verify that the regridding is equal: this is done by comparing the output from CDO to the output obtained by SMM. \n",
-    "The files to be checked are above. This is the same as the thing done in the tests."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "tas-healpix2.nc: remapnn via Dataset -> True\n",
-      "tas-healpix2.nc: remapnn via DataArray -> True\n",
-      "tas-healpix2.nc: remapcon via Dataset -> True\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "Ignoring index file '/home/paolo/smmregrid/tests/data/lsm-ifs.grb.923a8.idx' older than GRIB file\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "tas-healpix2.nc: remapcon via DataArray -> True\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "Ignoring index file '/home/paolo/smmregrid/tests/data/lsm-ifs.grb.923a8.idx' older than GRIB file\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "lsm-ifs.grb: remapnn via Dataset -> True\n",
-      "lsm-ifs.grb: remapnn via DataArray -> True\n",
-      "onlytos-ipsl.nc: remapnn via Dataset -> True\n",
-      "onlytos-ipsl.nc: remapnn via DataArray -> True\n",
-      "onlytos-ipsl.nc: remapcon via Dataset -> True\n",
-      "onlytos-ipsl.nc: remapcon via DataArray -> True\n",
-      "onlytos-ipsl.nc: remapbil via Dataset -> True\n",
-      "onlytos-ipsl.nc: remapbil via DataArray -> True\n",
-      "tas-ecearth.nc: remapnn via Dataset -> True\n",
-      "tas-ecearth.nc: remapnn via DataArray -> True\n",
-      "tas-ecearth.nc: remapcon via Dataset -> True\n",
-      "tas-ecearth.nc: remapcon via DataArray -> True\n",
-      "tas-ecearth.nc: remapbil via Dataset -> True\n",
-      "tas-ecearth.nc: remapbil via DataArray -> True\n",
-      "2t-era5.nc: remapnn via Dataset -> True\n",
-      "2t-era5.nc: remapnn via DataArray -> True\n",
-      "2t-era5.nc: remapcon via Dataset -> True\n",
-      "2t-era5.nc: remapcon via DataArray -> True\n",
-      "2t-era5.nc: remapbil via Dataset -> True\n",
-      "2t-era5.nc: remapbil via DataArray -> True\n",
-      "tos-fesom.nc: remapnn via Dataset -> True\n",
-      "tos-fesom.nc: remapnn via DataArray -> True\n",
-      "tos-fesom.nc: remapcon via Dataset -> True\n",
-      "tos-fesom.nc: remapcon via DataArray -> True\n",
-      "ua-ecearth.nc: remapnn via Dataset -> True\n",
-      "ua-ecearth.nc: remapnn via DataArray -> True\n",
-      "ua-ecearth.nc: remapcon via Dataset -> True\n",
-      "ua-ecearth.nc: remapcon via DataArray -> True\n",
-      "ua-ecearth.nc: remapbil via Dataset -> True\n",
-      "ua-ecearth.nc: remapbil via DataArray -> True\n",
-      "mix-cesm.nc: remapnn via Dataset -> True\n",
-      "mix-cesm.nc: remapcon via Dataset -> True\n",
-      "mix-cesm.nc: remapbil via Dataset -> True\n"
-     ]
-    }
-   ],
-   "source": [
-    "\n",
-    "for filein in filelist: \n",
-    "    for method in base[filein]['methods']:\n",
-    "        for access in base[filein]['accesses']: \n",
-    "            cc = check_cdo_regrid(os.path.join(indir,filein), tfile, method = method, access = access)#, extra = base[filein].get('extra'))\n",
-    "            print(filein + ': remap' + method + ' via ' + access + ' -> ' + str(cc))\n"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "From now we use only cases where we can use remapcon, i.e. we remove gaussiam reduced"
-   ]
-  },
   {
    "attachments": {},
    "cell_type": "markdown",
@@ -181,16 +78,89 @@
    "metadata": {},
    "outputs": [
     {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "caught signal <cdo.Cdo.__getattr__.<locals>.Operator object at 0x7fe995fca620> 2 <frame at 0x7fe8c829c040, file '/work/users/paolo/miniconda3/envs/smmregrid/lib/python3.10/selectors.py', line 416, code select>\n"
-     ]
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>CDO (Weight+Remap)</th>\n",
+       "      <th>SMM (Weight+Remap)</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>tas-healpix2.nc</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.801709</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>onlytos-ipsl.nc</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.870045</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>tas-ecearth.nc</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.753496</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2t-era5.nc</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.626265</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>tos-fesom.nc</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.924330</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>ua-ecearth.nc</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.719938</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>mix-cesm.nc</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.665692</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                 CDO (Weight+Remap)  SMM (Weight+Remap)\n",
+       "tas-healpix2.nc                 1.0            0.801709\n",
+       "onlytos-ipsl.nc                 1.0            0.870045\n",
+       "tas-ecearth.nc                  1.0            0.753496\n",
+       "2t-era5.nc                      1.0            0.626265\n",
+       "tos-fesom.nc                    1.0            0.924330\n",
+       "ua-ecearth.nc                   1.0            0.719938\n",
+       "mix-cesm.nc                     1.0            0.665692"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
     "# nrepetition for the check\n",
-    "nr = 5\n",
+    "nr = 10\n",
     "\n",
     "if 'lsm-ifs.grb' in filelist:\n",
     "    base.pop('lsm-ifs.grb')\n",
@@ -201,7 +171,6 @@
     "    xfield = xr.open_mfdataset(ifile)\n",
     "    wfield = cdo_generate_weights(ifile, tfile, method = 'con')\n",
     "    interpolator = Regridder(weights=wfield)\n",
-    "    var = list(xfield.data_vars)[-1]\n",
     "    rfield = interpolator.regrid(xfield)\n",
     "    return(rfield)\n",
     "\n",
@@ -226,27 +195,14 @@
    "source": [
     "# Remapping (with weights available)\n",
     "\n",
-    "This is the real goal of smmregrid. Here we test the computation of the remap when the weights are pre-computed, still using with conservative remapping. Considering that SMM does not have to write anything to disk, it is a few times faster. Running with Dataset implies a bit of overhead (20%). Masks have been integrated and create a small overhead when needed. Of course, loading the files implies a considerable slowdown."
+    "This is the real goal of smmregrid. Here we test the computation of the remap when the weights are pre-computed, still using with conservative remapping. Considering that SMM does not have to write anything to disk, it is a few times faster. Running with Dataset implies a bit of overhead (20%). Masks have been integrated and create a small overhead when needed. Of course, loading the files into memory implies a considerable slowdown."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "tas-healpix2.nc\n",
-      "onlytos-ipsl.nc\n",
-      "tas-ecearth.nc\n",
-      "2t-era5.nc\n",
-      "tos-fesom.nc\n",
-      "ua-ecearth.nc\n",
-      "mix-cesm.nc\n"
-     ]
-    },
     {
      "data": {
       "text/html": [
@@ -285,84 +241,84 @@
        "      <td>1</td>\n",
        "      <td>(12, 12288)</td>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.822694</td>\n",
-       "      <td>0.255584</td>\n",
-       "      <td>0.173050</td>\n",
-       "      <td>0.036148</td>\n",
-       "      <td>0.269337</td>\n",
-       "      <td>0.244628</td>\n",
+       "      <td>0.965832</td>\n",
+       "      <td>0.128650</td>\n",
+       "      <td>0.103219</td>\n",
+       "      <td>0.021833</td>\n",
+       "      <td>0.115453</td>\n",
+       "      <td>0.122626</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>onlytos-ipsl.nc</th>\n",
        "      <td>1</td>\n",
        "      <td>(12, 332, 362)</td>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.972039</td>\n",
-       "      <td>0.400861</td>\n",
-       "      <td>0.323148</td>\n",
-       "      <td>0.080454</td>\n",
-       "      <td>0.440603</td>\n",
-       "      <td>0.397795</td>\n",
+       "      <td>0.980560</td>\n",
+       "      <td>0.133137</td>\n",
+       "      <td>0.133919</td>\n",
+       "      <td>0.041179</td>\n",
+       "      <td>0.144053</td>\n",
+       "      <td>0.142921</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>tas-ecearth.nc</th>\n",
        "      <td>1</td>\n",
        "      <td>(12, 256, 512)</td>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.941557</td>\n",
-       "      <td>0.370303</td>\n",
-       "      <td>0.304809</td>\n",
-       "      <td>0.078599</td>\n",
-       "      <td>0.394931</td>\n",
-       "      <td>0.379545</td>\n",
+       "      <td>0.993987</td>\n",
+       "      <td>0.144356</td>\n",
+       "      <td>0.145000</td>\n",
+       "      <td>0.037764</td>\n",
+       "      <td>0.156681</td>\n",
+       "      <td>0.155994</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2t-era5.nc</th>\n",
        "      <td>1</td>\n",
        "      <td>(12, 73, 144)</td>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.982425</td>\n",
-       "      <td>0.228457</td>\n",
-       "      <td>0.160097</td>\n",
-       "      <td>0.041386</td>\n",
-       "      <td>0.259716</td>\n",
-       "      <td>0.258798</td>\n",
+       "      <td>1.102008</td>\n",
+       "      <td>0.087067</td>\n",
+       "      <td>0.079244</td>\n",
+       "      <td>0.022769</td>\n",
+       "      <td>0.105055</td>\n",
+       "      <td>0.092556</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>tos-fesom.nc</th>\n",
        "      <td>1</td>\n",
        "      <td>(12, 126859)</td>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.982059</td>\n",
-       "      <td>0.320393</td>\n",
-       "      <td>0.309535</td>\n",
-       "      <td>0.044901</td>\n",
-       "      <td>0.363529</td>\n",
-       "      <td>0.343803</td>\n",
+       "      <td>0.981255</td>\n",
+       "      <td>0.141750</td>\n",
+       "      <td>0.135510</td>\n",
+       "      <td>0.027498</td>\n",
+       "      <td>0.144329</td>\n",
+       "      <td>0.143400</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>ua-ecearth.nc</th>\n",
        "      <td>1</td>\n",
        "      <td>(2, 19, 256, 512)</td>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.897041</td>\n",
-       "      <td>0.519503</td>\n",
-       "      <td>0.468576</td>\n",
-       "      <td>0.152225</td>\n",
-       "      <td>0.707746</td>\n",
-       "      <td>0.695166</td>\n",
+       "      <td>0.986048</td>\n",
+       "      <td>0.236490</td>\n",
+       "      <td>0.233315</td>\n",
+       "      <td>0.068496</td>\n",
+       "      <td>0.251854</td>\n",
+       "      <td>0.271568</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>mix-cesm.nc</th>\n",
        "      <td>4</td>\n",
        "      <td>(12, 192, 288)</td>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.858966</td>\n",
-       "      <td>0.636257</td>\n",
-       "      <td>0.172329</td>\n",
-       "      <td>0.042839</td>\n",
-       "      <td>0.883419</td>\n",
-       "      <td>0.226711</td>\n",
+       "      <td>0.921705</td>\n",
+       "      <td>0.223586</td>\n",
+       "      <td>0.077339</td>\n",
+       "      <td>0.021130</td>\n",
+       "      <td>0.254226</td>\n",
+       "      <td>0.089087</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
@@ -370,49 +326,42 @@
       ],
       "text/plain": [
        "                 NVars           NRecords  CDO  CDO (NoLoad)  SMM (Dataset)  \\\n",
-       "tas-healpix2.nc      1        (12, 12288)  1.0      0.822694       0.255584   \n",
-       "onlytos-ipsl.nc      1     (12, 332, 362)  1.0      0.972039       0.400861   \n",
-       "tas-ecearth.nc       1     (12, 256, 512)  1.0      0.941557       0.370303   \n",
-       "2t-era5.nc           1      (12, 73, 144)  1.0      0.982425       0.228457   \n",
-       "tos-fesom.nc         1       (12, 126859)  1.0      0.982059       0.320393   \n",
-       "ua-ecearth.nc        1  (2, 19, 256, 512)  1.0      0.897041       0.519503   \n",
-       "mix-cesm.nc          4     (12, 192, 288)  1.0      0.858966       0.636257   \n",
+       "tas-healpix2.nc      1        (12, 12288)  1.0      0.965832       0.128650   \n",
+       "onlytos-ipsl.nc      1     (12, 332, 362)  1.0      0.980560       0.133137   \n",
+       "tas-ecearth.nc       1     (12, 256, 512)  1.0      0.993987       0.144356   \n",
+       "2t-era5.nc           1      (12, 73, 144)  1.0      1.102008       0.087067   \n",
+       "tos-fesom.nc         1       (12, 126859)  1.0      0.981255       0.141750   \n",
+       "ua-ecearth.nc        1  (2, 19, 256, 512)  1.0      0.986048       0.236490   \n",
+       "mix-cesm.nc          4     (12, 192, 288)  1.0      0.921705       0.223586   \n",
        "\n",
        "                 SMM (DataArray)  SMM (DataArray+NoLoad)  SMM (Dataset+Write)  \\\n",
-       "tas-healpix2.nc         0.173050                0.036148             0.269337   \n",
-       "onlytos-ipsl.nc         0.323148                0.080454             0.440603   \n",
-       "tas-ecearth.nc          0.304809                0.078599             0.394931   \n",
-       "2t-era5.nc              0.160097                0.041386             0.259716   \n",
-       "tos-fesom.nc            0.309535                0.044901             0.363529   \n",
-       "ua-ecearth.nc           0.468576                0.152225             0.707746   \n",
-       "mix-cesm.nc             0.172329                0.042839             0.883419   \n",
+       "tas-healpix2.nc         0.103219                0.021833             0.115453   \n",
+       "onlytos-ipsl.nc         0.133919                0.041179             0.144053   \n",
+       "tas-ecearth.nc          0.145000                0.037764             0.156681   \n",
+       "2t-era5.nc              0.079244                0.022769             0.105055   \n",
+       "tos-fesom.nc            0.135510                0.027498             0.144329   \n",
+       "ua-ecearth.nc           0.233315                0.068496             0.251854   \n",
+       "mix-cesm.nc             0.077339                0.021130             0.254226   \n",
        "\n",
        "                 SMM (DataArray+Write)  \n",
-       "tas-healpix2.nc               0.244628  \n",
-       "onlytos-ipsl.nc               0.397795  \n",
-       "tas-ecearth.nc                0.379545  \n",
-       "2t-era5.nc                    0.258798  \n",
-       "tos-fesom.nc                  0.343803  \n",
-       "ua-ecearth.nc                 0.695166  \n",
-       "mix-cesm.nc                   0.226711  "
+       "tas-healpix2.nc               0.122626  \n",
+       "onlytos-ipsl.nc               0.142921  \n",
+       "tas-ecearth.nc                0.155994  \n",
+       "2t-era5.nc                    0.092556  \n",
+       "tos-fesom.nc                  0.143400  \n",
+       "ua-ecearth.nc                 0.271568  \n",
+       "mix-cesm.nc                   0.089087  "
       ]
      },
-     "execution_count": 4,
+     "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "if 'lsm-ifs.grb' in filelist:\n",
-    "    base.pop('lsm-ifs.grb')\n",
-    "\n",
     "data =[]\n",
     "for filein in base.keys(): \n",
-    "    print(filein)\n",
-    "    if filein in ['era5-mon.nc', 'era5-day.nc']: \n",
-    "        nr = 2\n",
-    "    else :\n",
-    "        nr = 10\n",
+    "    nr = 10\n",
     "\n",
     "    # CDO\n",
     "    wfile = cdo.gencon(tfile, input = os.path.join(indir,filein))\n",
@@ -438,6 +387,7 @@
     "    arr =       timeit.timeit(lambda: interpolator.regrid(xfield[myvar[0]]).load(), number = nr)\n",
     "    arrnoload = timeit.timeit(lambda: interpolator.regrid(xfield[myvar[0]]), number = nr)\n",
     "    #arrnomask = timeit.timeit(lambda: interpolator.regrid(xfield[myvar[0]], masked = False).load(), number = nr)\n",
+    "    \n",
     "    setwrite =  timeit.timeit(lambda: interpolator.regrid(xfield).to_netcdf('test.nc'), number = nr)\n",
     "    if os.path.isfile('test.nc'):\n",
     "        os.remove('test.nc')\n",
@@ -455,39 +405,6 @@
     "final\n"
    ]
   },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Print to markdown so that we can copy paste it on the portal"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "|                 |   NVars | NRecords          |   CDO |   CDO (NoLoad) |   SMM (Dataset) |   SMM (DataArray) |   SMM (DataArray+NoLoad) |   SMM (Dataset+Write) |   SMM (DataArray+Write) |\n",
-      "|:----------------|--------:|:------------------|------:|---------------:|----------------:|------------------:|-------------------------:|----------------------:|------------------------:|\n",
-      "| tas-healpix2.nc |       1 | (12, 12288)       |     1 |       0.822694 |        0.255584 |          0.17305  |                0.0361485 |              0.269337 |                0.244628 |\n",
-      "| onlytos-ipsl.nc |       1 | (12, 332, 362)    |     1 |       0.972039 |        0.400861 |          0.323148 |                0.080454  |              0.440603 |                0.397795 |\n",
-      "| tas-ecearth.nc  |       1 | (12, 256, 512)    |     1 |       0.941557 |        0.370303 |          0.304809 |                0.0785989 |              0.394931 |                0.379545 |\n",
-      "| 2t-era5.nc      |       1 | (12, 73, 144)     |     1 |       0.982425 |        0.228457 |          0.160097 |                0.0413857 |              0.259716 |                0.258798 |\n",
-      "| tos-fesom.nc    |       1 | (12, 126859)      |     1 |       0.982059 |        0.320393 |          0.309535 |                0.0449009 |              0.363529 |                0.343803 |\n",
-      "| ua-ecearth.nc   |       1 | (2, 19, 256, 512) |     1 |       0.897041 |        0.519503 |          0.468576 |                0.152225  |              0.707746 |                0.695166 |\n",
-      "| mix-cesm.nc     |       4 | (12, 192, 288)    |     1 |       0.858966 |        0.636257 |          0.172329 |                0.0428392 |              0.883419 |                0.226711 |\n"
-     ]
-    }
-   ],
-   "source": [
-    "print(final.to_markdown())"
-   ]
-  },
   {
    "attachments": {},
    "cell_type": "markdown",
@@ -495,12 +412,12 @@
    "source": [
     "## Weight generation\n",
     "\n",
-    "Test the different weights generation possibilities with CDO, tested with conservative remapping: the climtas code is way more efficient if files are already on the disk, since the call to CDO has to be done from file. CDO bindings have a minimum overhead to be considered"
+    "As a final note, test the different weights generation possibilities with CDO, tested with conservative remapping: the climtas code is way more efficient if files are already on the disk, since the call to CDO has to be done from file. CDO bindings have a minimum overhead to be considered"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
@@ -531,40 +448,46 @@
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
+       "      <th>tas-healpix2.nc</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.694908</td>\n",
+       "      <td>0.723904</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
        "      <th>onlytos-ipsl.nc</th>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.938895</td>\n",
-       "      <td>1.025188</td>\n",
+       "      <td>0.837255</td>\n",
+       "      <td>0.879399</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>tas-ecearth.nc</th>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.882129</td>\n",
-       "      <td>1.164144</td>\n",
+       "      <td>0.709471</td>\n",
+       "      <td>0.892979</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2t-era5.nc</th>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.747073</td>\n",
-       "      <td>0.852617</td>\n",
+       "      <td>0.602283</td>\n",
+       "      <td>0.658382</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>tos-fesom.nc</th>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.987944</td>\n",
-       "      <td>1.377073</td>\n",
+       "      <td>0.902827</td>\n",
+       "      <td>1.253558</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>ua-ecearth.nc</th>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.891298</td>\n",
-       "      <td>1.118633</td>\n",
+       "      <td>0.702836</td>\n",
+       "      <td>0.805423</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>mix-cesm.nc</th>\n",
        "      <td>1.0</td>\n",
-       "      <td>0.796245</td>\n",
-       "      <td>1.012632</td>\n",
+       "      <td>0.646993</td>\n",
+       "      <td>0.737864</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
@@ -572,42 +495,32 @@
       ],
       "text/plain": [
        "                 CDO bindings  CDO subprocess (from file)  \\\n",
-       "onlytos-ipsl.nc           1.0                    0.938895   \n",
-       "tas-ecearth.nc            1.0                    0.882129   \n",
-       "2t-era5.nc                1.0                    0.747073   \n",
-       "tos-fesom.nc              1.0                    0.987944   \n",
-       "ua-ecearth.nc             1.0                    0.891298   \n",
-       "mix-cesm.nc               1.0                    0.796245   \n",
+       "tas-healpix2.nc           1.0                    0.694908   \n",
+       "onlytos-ipsl.nc           1.0                    0.837255   \n",
+       "tas-ecearth.nc            1.0                    0.709471   \n",
+       "2t-era5.nc                1.0                    0.602283   \n",
+       "tos-fesom.nc              1.0                    0.902827   \n",
+       "ua-ecearth.nc             1.0                    0.702836   \n",
+       "mix-cesm.nc               1.0                    0.646993   \n",
        "\n",
        "                 CDO subprocess (from xarray)  \n",
-       "onlytos-ipsl.nc                      1.025188  \n",
-       "tas-ecearth.nc                       1.164144  \n",
-       "2t-era5.nc                           0.852617  \n",
-       "tos-fesom.nc                         1.377073  \n",
-       "ua-ecearth.nc                        1.118633  \n",
-       "mix-cesm.nc                          1.012632  "
+       "tas-healpix2.nc                      0.723904  \n",
+       "onlytos-ipsl.nc                      0.879399  \n",
+       "tas-ecearth.nc                       0.892979  \n",
+       "2t-era5.nc                           0.658382  \n",
+       "tos-fesom.nc                         1.253558  \n",
+       "ua-ecearth.nc                        0.805423  \n",
+       "mix-cesm.nc                          0.737864  "
       ]
      },
-     "execution_count": 3,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
-    },
-    {
-     "ename": "",
-     "evalue": "",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[1;31mnotebook controller is DISPOSED. \n",
-      "\u001b[1;31mPer ulteriori dettagli, visualizzare Jupyter <a href='command:jupyter.viewOutput'>log</a>."
-     ]
     }
    ],
    "source": [
     "# nrepetition for the check\n",
-    "nr = 5\n",
-    "\n",
-    "if 'lsm-ifs.grb' in filelist:\n",
-    "    base.pop('lsm-ifs.grb')\n",
+    "nr = 10\n",
     "\n",
     "# generate weights from file\n",
     "data = []\n",
@@ -650,7 +563,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.8"
+   "version": "3.12.7"
   },
   "orig_nbformat": 4,
   "vscode": {


Test patch:
diff --git a/.github/workflows/mambatest.yml b/.github/workflows/mambatest.yml
index b4f6fb8..901d23e 100644
--- a/.github/workflows/mambatest.yml
+++ b/.github/workflows/mambatest.yml
@@ -21,7 +21,8 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        #python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.12"]
     defaults:
       run:
         shell: bash -el {0}
@@ -51,35 +52,6 @@ jobs:
         flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
         # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
-    - name: Install pytest
-      run: |
-         python -m pip install pytest
-    - name: Run tests
-      run: |
-        python -m pytest -v
-        
-  coverage:
-    needs: build
-    runs-on: ubuntu-latest
-    defaults:
-      run:
-        shell: bash -el {0}
-    steps:
-    - uses: actions/checkout@v4
-    - name: provision-with-micromamba
-      uses: mamba-org/setup-micromamba@v1
-      with:
-        environment-file: environment.yml
-        environment-name: smmregrid
-        cache-downloads: true
-        cache-environment: true
-        condarc: |
-          channels:
-            - conda-forge
-    - name: Install smmregrid
-      run: |
-        # install package
-        python -m pip install --editable .
     - name: Install coverage tools
       run: |
         python -m pip install coverage
diff --git a/.github/workflows/pypitest.yml b/.github/workflows/pypitest.yml
index 9128dcb..bb429aa 100644
--- a/.github/workflows/pypitest.yml
+++ b/.github/workflows/pypitest.yml
@@ -29,7 +29,7 @@ jobs:
         environment-file: environment.yml
         environment-name: smmregrid
         cache-downloads: true
-        cache-environment: true
+        cache-environment: false
         condarc: |
           channels:
             - conda-forge
diff --git a/tests/areas_test.py b/tests/areas_test.py
new file mode 100644
index 0000000..ac8324d
--- /dev/null
+++ b/tests/areas_test.py
@@ -0,0 +1,52 @@
+"""Set of basic tests for smmregrid"""
+
+import os
+import pytest
+from smmregrid import CdoGenerate
+
+
+INDIR = 'tests/data'
+tfile = os.path.join(INDIR, 'r360x180.nc')
+EARTH_SURFACE = 5.101*1e8 #km2
+OCEAN_SURFACE = 3.6*1e8
+RELTOL = 0.02 # %2 percent margin
+TOLERANCE = EARTH_SURFACE * RELTOL
+
+@pytest.mark.parametrize("filename,shape,kind", [
+    ('2t-era5.nc', (73, 144), "global"),
+    ('onlytos-ipsl.nc', (332, 362), "global"),
+    ('tos-fesom.nc', (126859,), "ocean"),
+    ('tas-healpix2.nc', (12288,), "global"),
+    ('so3d-nemo.nc', (292,362), "global")
+])
+def test_basic_areas_source(filename, shape, kind):
+    """Basic test for area generation from source file"""
+    generator = CdoGenerate(os.path.join(INDIR, filename), tfile, loglevel='debug')
+    area = generator.areas()
+    compare_surface = EARTH_SURFACE if kind == "global" else OCEAN_SURFACE
+    earth_surface = area.cell_area.values.sum()/1e6
+    assert area.cell_area.shape == shape
+    assert earth_surface == pytest.approx(compare_surface, abs=TOLERANCE)
+
+@pytest.mark.parametrize("filename,shape", [
+    (os.path.join(INDIR, 'r360x180.nc'), (180, 360)),
+    ('hp32', (12288,)),
+    ('r180x90', (90, 180,)),
+    (os.path.join(INDIR, 'tas-healpix2.nc'), (12288,))
+])
+def test_basic_areas_target(filename, shape):
+    """Basic test for area generation from source file"""
+    infile = os.path.join(INDIR, 'r360x180.nc')
+    generator = CdoGenerate(infile, filename, loglevel='debug')
+    area = generator.areas(target=True)
+    earth_surface = area.cell_area.values.sum()/1e6
+    assert area.cell_area.shape == shape
+    assert earth_surface == pytest.approx(EARTH_SURFACE, abs=TOLERANCE)
+
+def test_nosource_areas_target():
+    """Basic test for area generation from source file"""
+    generator = CdoGenerate(source_grid=None, target_grid='r180x91', loglevel='debug')
+    area = generator.areas(target=True)
+    earth_surface = area.cell_area.values.sum()/1e6
+    assert area.cell_area.shape == (91, 180)
+    assert earth_surface == pytest.approx(EARTH_SURFACE, abs=TOLERANCE)
\ No newline at end of file
diff --git a/tests/basic_test.py b/tests/basic_test.py
index d7c7a61..61010bd 100644
--- a/tests/basic_test.py
+++ b/tests/basic_test.py
@@ -1,32 +1,64 @@
-"""Set of basic tests for smmregrid"""
+"""Set of basic tests for smmregrid weights"""
 
 import os
 import pytest
 import xarray
 import numpy
-from smmregrid.checker import cdo_generate_weights, Regridder
+from smmregrid import Regridder, CdoGenerate, cdo_generate_weights
+
 
 INDIR = 'tests/data'
 tfile = os.path.join(INDIR, 'r360x180.nc')
-rfile = os.path.join(INDIR, 'regional.nc')
 
 @pytest.mark.parametrize("method", ['con', 'nn', 'bil'])
 def test_healpix_extra(method):
     """Test for healpix with cdo_extra and cdo_options"""
-    wfield = cdo_generate_weights(os.path.join(INDIR, 'healpix_0.nc'), tfile,
-                                  method = method, cdo_extra = '-setgrid,hp1_nested',
-                                  cdo_options='--force', loglevel='debug')
+    if method == 'con':
+        options = ['--force']
+    else:
+        options = ['--force', '-f', 'nc']
+    wfield = CdoGenerate(os.path.join(INDIR, 'healpix_0.nc'), tfile,
+                                  cdo_extra = '-setgrid,hp1_nested',
+                                  cdo_options=options,
+                                  loglevel='debug').weights(method = method)
     interpolator = Regridder(weights=wfield, loglevel='debug')
     xfield = xarray.open_mfdataset(os.path.join(INDIR, 'healpix_0.nc'))
     rfield = interpolator.regrid(xfield)
     assert rfield['tas'].shape == (2, 180, 360)
 
-@pytest.mark.parametrize("method", ['con', 'nn', 'bil'])
+@pytest.mark.parametrize("method", ['con', 'nn', 'bic'])
 def test_nan_preserve(method):
-    """Test to verify that NaN are preserved"""
+    """Test to verify that NaN are preserved with deprecated cdo_generate_weights"""
     xfield = xarray.open_mfdataset(os.path.join(INDIR, 'tas-ecearth.nc'))
     xfield['tas'][1,:,:] = numpy.nan
-    wfield = cdo_generate_weights(xfield, tfile, method = method, loglevel='debug')
-    interpolator = Regridder(weights=wfield, loglevel='debug')
+    wfield = cdo_generate_weights(xfield, tfile, loglevel='debug', method = method)
+    interpolator = Regridder(weights=wfield, space_dims='pippo', loglevel='debug')
     rfield = interpolator.regrid(xfield)
     assert numpy.isnan(rfield['tas'][1,:,:]).all().compute()
+
+@pytest.mark.parametrize("method", ['nn'])
+def test_datarray(method):
+    """"Minimal test to verify regridding from DataArray"""
+    xfield = xarray.open_mfdataset(os.path.join(INDIR, 'tas-ecearth.nc'))
+    interpolator = Regridder(source_grid=xfield['tas'], target_grid=tfile,
+                             loglevel='debug', method = method)
+    interp = interpolator.regrid(source_data=xfield)
+    assert interp['tas'].shape == (12, 180, 360)
+
+@pytest.mark.parametrize("method", ['dis', 'con'])
+def test_horizontal_dims(method):
+    """"Minimal test to verify regridding with horizontal_dims coords"""
+    xfield = xarray.open_mfdataset(os.path.join(INDIR, 'tas-ecearth.nc'))
+    interpolator = Regridder(source_grid=xfield, target_grid=tfile, loglevel='debug',
+                             method = method, horizontal_dims=['lon', 'lat'])
+    interp = interpolator.regrid(source_data=xfield.isel(time=0))
+    assert interp['tas'].shape == (180, 360)
+
+@pytest.mark.parametrize("method", ['bil', 'con'])
+def test_toward_healpix(method):
+    """"Testing toward healpix string-defined-grid grid"""
+    xfield = xarray.open_mfdataset(os.path.join(INDIR, 'tas-ecearth.nc'))
+    interpolator = Regridder(source_grid=xfield, target_grid='hp16_nested', loglevel='debug',
+                             method = method, cdo_options='--force')
+    interp = interpolator.regrid(source_data=xfield)
+    assert interp['tas'].shape == (12, 3072)
diff --git a/tests/data/ua-so_mix_ecearth.nc b/tests/data/ua-so_mix_ecearth.nc
new file mode 100644
index 0000000..9e9f589
Binary files /dev/null and b/tests/data/ua-so_mix_ecearth.nc differ
diff --git a/tests/gridinspector_test.py b/tests/gridinspector_test.py
new file mode 100644
index 0000000..3969db8
--- /dev/null
+++ b/tests/gridinspector_test.py
@@ -0,0 +1,25 @@
+"""Set of basic tests for smmregrid gridinspector"""
+
+import os
+import pytest
+import xarray
+from smmregrid import GridInspector
+
+
+INDIR = 'tests/data'
+
+@pytest.mark.parametrize("file,ngrids,firstdims, variables", [
+    ("2t-era5.nc", 1, ["lon", "lat"], ['2t']),
+    ("mix-cesm.nc", 1, ["lon", "lat"], ['hfss', 'hfls', 'rlds', 'rlus']),
+    ("so3d-nemo.nc", 1, ["i", "j", "lev"], ['so']),
+    ("ua-so_mix_ecearth.nc", 2, ["lon", "lat"], ['ua']),
+    ("temp3d-fesom.nc", 1, ['nod2', 'nz1'], ['temp']),
+    ])
+def test_basic_gridinspector(file, ngrids, firstdims, variables):
+    """test for GridInspector"""
+    filename = os.path.join('tests/data', file)
+    xfield = xarray.open_mfdataset(filename)
+    grids = GridInspector(xfield, loglevel='debug').get_grid_info()
+    assert len(grids) == ngrids
+    assert set(grids[0].dims) == set(firstdims)
+    assert set(grids[0].variables.keys()) == set(variables)
\ No newline at end of file
diff --git a/tests/gridtype_test.py b/tests/gridtype_test.py
new file mode 100644
index 0000000..d58dc75
--- /dev/null
+++ b/tests/gridtype_test.py
@@ -0,0 +1,41 @@
+"""Set of basic tests for smmregrid weights"""
+
+import pytest
+from smmregrid import GridType
+
+@pytest.mark.parametrize("definition,dims,other", [
+    (['lon','lat'], ['lon', 'lat'], []),
+    (['lon','lat', 'lev', 'time'], ['lon', 'lat', 'lev'], []),
+    (['i', 'k'], ['i'], ['k']),
+    (['pix', 'time', 'papera'], ['pix'], ['papera']),
+    (['time', 'papera'], [], ['papera'])
+
+    ])
+def test_gridtype(definition, dims, other):
+    """Basic gridtype investigation"""
+    grid = GridType(dims=definition)
+    assert set(grid.dims) == set(dims)
+    assert set(grid.other_dims) == set(other)
+
+def test_difference():
+    """Test equality for gridtypes"""
+    grid1 = GridType(dims=['lon', 'lat'])
+    grid2 = GridType(dims=['i', 'j'])
+    grid3 = GridType(dims=['lon', 'lat', 'time', 'plev'])
+    assert grid1 != grid2
+    assert grid1 == grid3
+
+def test_multiple_vertical():
+    """Test multiple vertical"""
+    with pytest.raises(ValueError):
+        GridType(dims=['lon', 'lat', 'lev', 'nz1'])
+
+def test_gridtype_extradims():
+    """Test for extradimensions"""
+    #with pytest.raises(ValueError):
+    #    GridType(dims=["lon", "lat", "ciccio"],
+    #             extra_dims = {'horizontal': "nonnapaper"})
+
+    grid = GridType(dims=["lon", "lat", "ciccio"],
+                    extra_dims = {'vertical': ["ciccio"]})
+    assert grid.vertical_dim == "ciccio"
diff --git a/tests/identity_test.py b/tests/identity2d_test.py
similarity index 64%
rename from tests/identity_test.py
rename to tests/identity2d_test.py
index 3e13541..783adad 100644
--- a/tests/identity_test.py
+++ b/tests/identity2d_test.py
@@ -2,25 +2,13 @@
 
 import os
 import pytest
+import xarray as xr
 from smmregrid.checker import check_cdo_regrid
 
 INDIR = 'tests/data'
 tfile = os.path.join(INDIR, 'r360x180.nc')
 rfile = os.path.join(INDIR, 'regional.nc')
 
-# test for pressure levels on gaussian grid (3D)
-@pytest.mark.parametrize("method", ['con', 'nn'])
-def test_nemo_3d(method):
-    fff = check_cdo_regrid(os.path.join(INDIR, 'so3d-nemo.nc'), tfile,
-                           remap_method=method, init_method='grids')
-    assert fff is True
-
-# test for pressure levels on gaussian grid (3D)
-@pytest.mark.parametrize("method", ['con', 'nn'])
-def test_fesom_3d(method):
-    fff = check_cdo_regrid(os.path.join(INDIR, 'temp3d-fesom.nc'), tfile,
-                           remap_method=method, init_method='grids')
-    assert fff is True
 
 # test for gaussian reduced grid (only nn)
 @pytest.mark.parametrize("method", ['nn', 'con'])
@@ -37,21 +25,22 @@ def test_gaussian_reduced(method):
     assert fff is True
 
 # test for gaussian grids as EC-Earth cmor
-@pytest.mark.parametrize("method", ['bil', 'con', 'nn'])
+@pytest.mark.parametrize("method", ['con', 'nn'])
 def test_gaussian_regular(method):
-    fff = check_cdo_regrid(os.path.join(INDIR, 'tas-ecearth.nc'), tfile,
+    xfield = xr.open_dataset(os.path.join(INDIR, 'tas-ecearth.nc'))
+    fff = check_cdo_regrid(xfield, tfile,
                            remap_method=method)
     assert fff is True
 
 # test for gaussian grids as EC-Earth cmor
-@pytest.mark.parametrize("method", ['bil', 'con', 'nn'])
+@pytest.mark.parametrize("method", ['bil', 'con'])
 def test_gaussian_regular_regional(method):
     fff = check_cdo_regrid(os.path.join(INDIR, 'tas-ecearth.nc'), rfile,
                            remap_method=method)
     assert fff is True
 
 # test for lonlt grids, init by weights
-@pytest.mark.parametrize("method", ['bil', 'con', 'nn'])
+@pytest.mark.parametrize("method", [ 'con', 'bil'])
 def test_lonlat(method):
     fff = check_cdo_regrid(os.path.join(INDIR, '2t-era5.nc'), tfile,
                            remap_method=method, init_method='weights')
@@ -66,21 +55,14 @@ def test_unstructured(method):
 
 # test for curvilinear grid
 @pytest.mark.parametrize("method", ['con', 'nn', 'bil'])
-def test_curivilinear(method):
+def test_curvilinear(method):
     fff = check_cdo_regrid(os.path.join(INDIR, 'onlytos-ipsl.nc'), tfile,
                            remap_method=method, init_method='grids')
     assert fff is True
 
-# test for pressure levels on gaussian grid (2D, level-by-level), init by weights
-@pytest.mark.parametrize("method", ['con', 'nn', 'bil'])
-def test_levbylev_plev_gaussian(method):
-    fff = check_cdo_regrid(os.path.join(INDIR, 'ua-ecearth.nc'), tfile,
-                           remap_method=method, init_method='weights',
-                           vert_coord="plev")
-    assert fff is True
 
 # test for pressure levels on gaussian grid with info logging (3D)
-@pytest.mark.parametrize("method", ['con', 'nn', 'bil'])
+@pytest.mark.parametrize("method", ['con', 'bil'])
 def test_full_plev_gaussian(method):
     fff = check_cdo_regrid(os.path.join(INDIR, 'ua-ecearth.nc'), tfile,
                            remap_method=method, init_method='grids')
diff --git a/tests/identity3d_test.py b/tests/identity3d_test.py
new file mode 100644
index 0000000..0361d3b
--- /dev/null
+++ b/tests/identity3d_test.py
@@ -0,0 +1,36 @@
+"""Set of basic tests for smmregrid"""
+
+import os
+import pytest
+import xarray as xr
+from smmregrid.checker import check_cdo_regrid
+
+INDIR = 'tests/data'
+tfile = os.path.join(INDIR, 'r360x180.nc')
+rfile = os.path.join(INDIR, 'regional.nc')
+
+# test for NEMO 3d grid
+@pytest.mark.parametrize("method", ['con', 'nn'])
+def test_nemo_3d(method):
+    if method == 'nn':
+        xfield = os.path.join(INDIR, 'so3d-nemo.nc')
+    else:
+        xfield = xr.open_dataset(os.path.join(INDIR, 'so3d-nemo.nc'))
+    fff = check_cdo_regrid(xfield, tfile,
+                           remap_method=method, init_method='grids')
+    assert fff is True
+
+# test for for FESOM 3d grid
+@pytest.mark.parametrize("method", ['con'])
+def test_fesom_3d(method):
+    fff = check_cdo_regrid(os.path.join(INDIR, 'temp3d-fesom.nc'), tfile,
+                           remap_method=method, init_method='grids')
+    assert fff is True
+
+# test for pressure levels on gaussian grid (2D, level-by-level), init by weights
+@pytest.mark.parametrize("method", ['con'])
+def test_levbylev_plev_gaussian(method):
+    fff = check_cdo_regrid(os.path.join(INDIR, 'ua-ecearth.nc'), tfile,
+                           remap_method=method, init_method='weights',
+                           vertical_dim="plev")
+    assert fff is True
\ No newline at end of file
diff --git a/tests/multigrid_test.py b/tests/multigrid_test.py
new file mode 100644
index 0000000..f3dfe14
--- /dev/null
+++ b/tests/multigrid_test.py
@@ -0,0 +1,30 @@
+"""Tests for smmregrid to verify level selection"""
+
+import os
+import pytest
+import xarray as xr
+from smmregrid.checker import check_cdo_regrid, Regridder
+
+INDIR = 'tests/data'
+ifile = os.path.join(INDIR, 'ua-so_mix_ecearth.nc')
+tfile = os.path.join(INDIR, 'r360x180.nc')
+
+def test_multi_grid_inspection():
+    """test and assert with CDO multigrids"""
+
+    regrid = Regridder(ifile, 'r180x90', loglevel='debug')
+    data = xr.open_dataset(ifile)
+    regridded = regrid.regrid(data)
+    assert len(regrid.grids) == 2
+    assert sorted(vars(regrid.grids[1])['dims']) == sorted(['lev', 'j', 'i'])
+    assert sorted(vars(regrid.grids[0])['dims']) == sorted(['plev', 'lat', 'lon'])
+    assert regridded['ua'].shape == (2, 3, 90, 180)
+    assert regridded['so'].shape == (2, 3, 90, 180)
+
+
+def test_multi_grid_cdo():
+    """test and assert with CDO multigrids"""
+
+    assert check_cdo_regrid(ifile, tfile, init_method='grids')
+    with pytest.raises(ValueError):
+        assert check_cdo_regrid(ifile, tfile, init_method='weights')
\ No newline at end of file
diff --git a/tests/util_test.py b/tests/util_test.py
new file mode 100644
index 0000000..40b840a
--- /dev/null
+++ b/tests/util_test.py
@@ -0,0 +1,53 @@
+"""Tests for util module"""
+
+import os
+import pytest
+import xarray
+from smmregrid.util import is_cdo_grid, check_gridfile
+
+
+@pytest.mark.parametrize("grid_str,expected", [
+    ("global_1.0", True),        # Global regular grid
+    ("dcw:US", True),            # Regional regular grid
+    ("zonal_2.5", True),         # Zonal latitudes
+    ("r360x180", True),          # Global regular NxM grid
+    ("lon=-75.0/lat=40.0", True),# One grid point
+    ("F64", True),               # Full regular Gaussian grid
+    ("n400", True),              # Full regular Gaussian grid
+    ("gme10", True),             # Global icosahedral-hexagonal GME grid
+    ("hp1024", True),            # HEALPix grid
+    ("hp32_ring", True),         # HEALPix grid
+    ("hpz4", True),              # HEALPix zoom grid
+    ("random_string", False),    # Invalid CDO grid
+    ("/path/to/file.nc", False)  # Invalid CDO grid (file path)
+])
+def test_is_cdo_grid(grid_str, expected):
+    """Simple test for cdo grid detection"""
+    assert is_cdo_grid(grid_str) == expected
+
+# Mock function to avoid actual file system checks for testing purposes
+def mock_exists(path):
+    return path == "/valid/file/path.nc"
+
+# Patch os.path.exists to use the mock function
+@pytest.fixture(autouse=True)
+def mock_os_path_exists(monkeypatch):
+    monkeypatch.setattr(os.path, "exists", mock_exists)
+
+@pytest.mark.parametrize("filename,expected,raises", [
+    (None, None, None),                                  # None input
+    (xarray.Dataset(), "xarray", None),                      # xarray.Dataset input
+    (xarray.DataArray(), "xarray", None),                    # xarray.DataArray input
+    ("global_1.0", "grid", None),                        # CDO grid string
+    ("r360x180", "grid", None),                          # CDO grid string
+    ("/valid/file/path.nc", "file", None),               # Valid file path
+    ("/invalid/file/path.nc", None, FileNotFoundError),  # Invalid file path
+    (123, None, TypeError),                              # Invalid type (not supported)
+])
+def test_check_gridfile(filename, expected, raises):
+    """simple test for check grid file"""
+    if raises:
+        with pytest.raises(raises):
+            check_gridfile(filename)
+    else:
+        assert check_gridfile(filename) == expected