Breaking changes: change spherical harmonic index nomenclature to 'nm…

…' and consistently use negative orders for Sine coefficients. Also improve documentation, and decide on numpy doctring style

docs/source/api.rst

@@ -1,6 +1,6 @@
-Python API Reference
-====================
-
+API Reference
+=============
+The calling 
 .. toctree::
    :maxdepth: 2
 

docs/source/conf.py

@@ -36,7 +36,7 @@
 html_theme = 'sphinx_rtd_theme'
 html_static_path = ['_static']
 
-
+napoleon_numpy_docstring = True
 
 nbsphinx_prolog = """
 Download this Jupyter notebook from `github <https://github.com/ITC-Water-Resources/shxarray/blob/main/docs/source/notebooks/{{ env.doc2path(env.docname, base=None) }}>`_

docs/source/installation.rst

@@ -1,4 +1,32 @@
 Installation and General Usage
 ==============================
+The latest **shxarray** is hosted on `pypi <https://pypi.org/project/shxarray/>`_ and can be installed through pip: 
+
+``pip install shxarray`` 
+
+Part of the module is written is `Cython <https://cython.readthedocs.io/en/latest/>`_, which means that a c compiler is needed to build the package. A binary wheel is currently not offered, but this may be offered in the future.
+
+Import and usage
+----------------
+For most operations, a simple import will expose the xarray extensions. For example:
+
+.. code-block:: python
+   
+   import shxarray
+   import xarray as xr
+
+   #Initialize a dataarray with zeros which has a dimension spanning degrees from nmin to nmax
+
+   nmax=20
+   nmin=2
+   dazeros=xr.DataArray.sh.ones(nmax=nmax,nmin=nmin)
+
+
+
+Development
+-----------
+Users interested in developing can install the latest version from `github <https://github.com/ITC-Water-Resources/shxarray/tree/main>`_. Cython is needed in case the binary extension is being developed, and users can consult the dedicated instructions on the github repository.
+
+Code can be supplied with `numpy docstrings <https://www.sphinx-doc.org/en/master/usage/extensions/example_numpy.html>`_ so they can be parsed into this documentation.
+
 
-The latest **shxarray** can be installed through ``pip install shxarray`` 

docs/source/introduction.rst

@@ -1,4 +1,26 @@
 Introduction
 ============
 
-The shxarray package is aimed to make using  spherical harmonic operations more accessible to a community which is used to xarray.
+The **shxarray** package is aimed to make using spherical harmonic operations more accessible to a community which is used to xarray.
+
+Putting degrees and orders in a MultiIndex
+------------------------------------------
+
+Spherical harmonic coefficients are often put in 2-dimensional arrays with degree and order spanning the 2 dimensions. This has the advantage that individual coefficients can be easily referenced as e.g. `cnm=mat[n,m]`. However, since only the upper triangle of those matrices are non-zero, the sparseness of these matrices is not made use of. When working with large datasets which also span other dimensions such as time or different levels, this will cause large segments of zeros.
+
+A `pandas.MultiIndex <https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.MultiIndex.html>`_ can facilitate the stacking of degrees and orders in a single dimension. On top of that, this multindex can then be used in `xarray` to work with spherical harmonics along a single coordinate. In **shxarray**, the spherical harmonic index is generally denoted as ``nm``, while when two spherical harmonic coordinates are needed alternative versions such as ``nm_`` are added. 
+
+Exposing spherical harmonic functionality to xarray
+---------------------------------------------------
+
+Some of the functionality needed for working with spherical harmonics need specialized access to the degree and order information. The aim of **shxarray** is to expose common functionality through `xarray accessors <https://docs.xarray.dev/en/stable/internals/extending-xarray.html>`_. This allows, besides familiar syntax for xarray users, also chaining of operations in a compact syntax.
+
+Delegate common operations to xarray
+------------------------------------
+
+In contrast to specialized spherical harmonic operations, many operations on spherical harmonic data can be delegated to xarray itself. Wherever possible, functionality and broadcasting features of xarray are made use for a consistent syntax.
+
+
+
+
+

docs/source/notebooks/visualize_filter.ipynb

@@ -141,15 +141,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 6,
    "id": "8bd71cd5-91d0-40d3-83f6-b015f04f501b",
    "metadata": {},
    "outputs": [],
    "source": [
     "# Note setting truncate to False below keeps coefficients below the minimum degree of the filter to their original values\n",
     "ddk4=daunit.sh.filter('DDK4',truncate=False)\n",
     "# assign to dataset variable but make sure they are sorted in the same order\n",
-    "dsfiltered['ddk4']=ddk4.sel(shi=dsfiltered.shi)\n"
+    "dsfiltered['ddk4']=ddk4.sel(nm=dsfiltered.nm)\n"
    ]
   },
   {
@@ -162,7 +162,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 7,
    "id": "725bb29f-9cb7-43ae-aad3-647e6052589e",
    "metadata": {},
    "outputs": [],
@@ -179,7 +179,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 8,
    "id": "bb124eb7-e8cc-48f5-8a3b-1675ae054cd5",
    "metadata": {},
    "outputs": [],
@@ -206,7 +206,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 9,
    "id": "c84fcfa3-8b54-4606-9be7-dba5ea6b8a47",
    "metadata": {},
    "outputs": [
@@ -216,7 +216,7 @@
        "Text(0.5, 1.0, 'DDK at lon,lat: (-60,-5)')"
       ]
      },
-     "execution_count": 8,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -280,7 +280,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.7"
   }
  },
  "nbformat": 4,

src/builtin_backend/Ynm.hpp

@@ -32,7 +32,7 @@ template <class ftype>
 class Ynm_cpp{
    public:
     Ynm_cpp(const int nmax);
-    Ynm_cpp(const size_t size, const int n [],const int m[ ],const int t[]);
+    Ynm_cpp(const size_t size, const int n [],const int m[ ]);
     Ynm_cpp(){};
     void set(const ftype lon, const ftype lat);
 
@@ -81,16 +81,15 @@ Ynm_cpp<ftype>::Ynm_cpp(int nmax):legnm(nmax),sz_(2*(legnm.idx(nmax,nmax)+1)-(nm
 }
 
 template<class ftype>
-Ynm_cpp<ftype>::Ynm_cpp(const size_t size, const int n [],const int m[ ],const int t[]):sz_(size),ynmdata_(size,0.0),mnidx_(sz_){
+Ynm_cpp<ftype>::Ynm_cpp(const size_t size, const int n [],const int m[ ]):sz_(size),ynmdata_(size,0.0),mnidx_(sz_){
 
     ///find nmax
     int nmax=-1;
     //create a temporary map used to figure out the correct indices of the desired output order
     //map<std::pair<int,int>,ssize_t> mninputidx;
     for (size_t i=0;i<size;++i){
 	nmax=(nmax<n[i])?n[i]:nmax;
-	int sgn=t[i]?-1:1;
-	mnidx_[i]={sgn*m[i],n[i],i};
+	mnidx_[i]={m[i],n[i],i};
     }
 
     legnm=Legendre_nm<ftype>(nmax);

src/builtin_backend/analysis.pyx

@@ -27,7 +27,7 @@ cdef class Analysis:
     def __cinit__(self, int nmax):
 
         #create a spherical harmonic index
-        self._dsobj=xr.Dataset(coords=SHindexBase.shi(nmax,0))
+        self._dsobj=xr.Dataset(coords=SHindexBase.nm(nmax,0))
 
     def __call__(self,dain:xr.DataArray):
         """Perform  spherical harmonic analysis on an input xarray DataArray object""" 
@@ -72,17 +72,16 @@ cdef class Analysis:
 
         cdef int auxsize=np.prod([val for ky,val in dain.sizes.items() if ky not in ["lon","lat"]])
 
-        cdef int shsize=len(self._dsobj.indexes['shi'])
+        cdef int shsize=len(self._dsobj.indexes['nm'])
         #memoryview to output data (sh dimension should vary quickest) 
         cdef double [:,:] outv=dout.data.reshape([auxsize,shsize])
         #This is the same a s a Fortran contiguous array with dimension shsize,auxsize, lada=shsize  
 
 
 
 
-        cdef int[::1] nv = self._dsobj.shi.n.data.astype(np.int32)
-        cdef int[::1] mv = self._dsobj.shi.m.data.astype(np.int32)
-        cdef int[::1] tv = self._dsobj.shi.t.data.astype(np.int32)
+        cdef int[::1] nv = self._dsobj.nm.n.data.astype(np.int32)
+        cdef int[::1] mv = self._dsobj.nm.m.data.astype(np.int32)
 
         cdef Ynm_cpp[double] ynm
         cdef int ilat,ilon          
@@ -126,7 +125,7 @@ cdef class Analysis:
         omp_init_lock(&lock)
         with nogil, parallel():
 
-            ynm=Ynm_cpp[double](shsize,&nv[0],&mv[0],&tv[0])
+            ynm=Ynm_cpp[double](shsize,&nv[0],&mv[0])
             for ilat in prange(nlat):
                 alpha=weight*cos(latv[ilat]*d2r)
                 for ilon in range(nlon):

src/builtin_backend/legendre.pxd

@@ -42,9 +42,9 @@ cdef extern from "Ynm.hpp":
     cdef cppclass Ynm_cpp[T] nogil:
         Ynm_cpp() except +
         Ynm_cpp(int nmax) except +
-        Ynm_cpp(cython.size_t size, const int n[],const int m[], const int t[]) except +
+        Ynm_cpp(cython.size_t size, const int n[],const int m[]) except +
         void set( T lon, T lat) nogil 
-        cython.ssize_t idx(int n,int m,int t)
+        cython.ssize_t idx(int n,int m)
         T& operator[](size_t i)
         int nmax()
         T* data()