Skip to content

Commit

Permalink
Breaking changes: change spherical harmonic index nomenclature to 'nm…
Browse files Browse the repository at this point in the history
…' and consistently use negative orders for Sine coefficients. Also improve documentation, and decide on numpy doctring style
  • Loading branch information
strawpants committed Feb 15, 2024
1 parent 1107a23 commit 994b43a
Show file tree
Hide file tree
Showing 24 changed files with 1,227 additions and 1,587 deletions.
6 changes: 3 additions & 3 deletions docs/source/api.rst
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Python API Reference
====================

API Reference
=============
The calling
.. toctree::
:maxdepth: 2

Expand Down
2 changes: 1 addition & 1 deletion docs/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -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) }}>`_
Expand Down
30 changes: 29 additions & 1 deletion docs/source/installation.rst
Original file line number Diff line number Diff line change
@@ -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``
24 changes: 23 additions & 1 deletion docs/source/introduction.rst
Original file line number Diff line number Diff line change
@@ -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.





100 changes: 54 additions & 46 deletions docs/source/notebooks/TerrestrialWaterStorage.ipynb

Large diffs are not rendered by default.

14 changes: 7 additions & 7 deletions docs/source/notebooks/visualize_filter.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
},
{
Expand All @@ -162,7 +162,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 7,
"id": "725bb29f-9cb7-43ae-aad3-647e6052589e",
"metadata": {},
"outputs": [],
Expand All @@ -179,7 +179,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 8,
"id": "bb124eb7-e8cc-48f5-8a3b-1675ae054cd5",
"metadata": {},
"outputs": [],
Expand All @@ -206,7 +206,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 9,
"id": "c84fcfa3-8b54-4606-9be7-dba5ea6b8a47",
"metadata": {},
"outputs": [
Expand All @@ -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"
},
Expand Down Expand Up @@ -280,7 +280,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.6"
"version": "3.11.7"
}
},
"nbformat": 4,
Expand Down
7 changes: 3 additions & 4 deletions src/builtin_backend/Ynm.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down Expand Up @@ -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);
Expand Down
11 changes: 5 additions & 6 deletions src/builtin_backend/analysis.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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"""
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down
4 changes: 2 additions & 2 deletions src/builtin_backend/legendre.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
Loading

0 comments on commit 994b43a

Please sign in to comment.