Monday, March 11, 2013

Spherical Harmonics transform and rotation in Python

Recently I needed a module that provides spherical harmonics transforms and tools in Python and I started wrapping the library shtools written by Mark Wieczorek using f2py from the numpy/scipy packages. Here is the result in case anyone is interested. If you are looking for a complete package, note that there is also another python module (pyspharm) that deals with spherical harmonics .

UPDATE: a new and better version of the python wrapper is now included in the SHTOOLS library which you can find here: pyshtools. Numerous examples are included that show how it can be used from Python. Other spherical harmonics libraries are: pyspharm and SHTns. In addition to the spherical harmonics transforms, SHTOOLS provides local multitaper spherical harmonics analysis, as well as spherical harmonics rotations, coupling matrices, etc ...

Little Examples

Once the package is compiled, a spherical harmonics transform is as easy as:
>>> import shtools
>>> grid = shtools.pymakegriddh(coeffs)
>>> coeffs = shtools.pyshexpanddh(grid,nlat=90)
a rotation by 3 Euler angles can be done with:
>>> angles = np.radians(np.array[90.,-30.,0.])) 
>>> coeffs = shtools.pyshrotaterealcoef(coeffs_trim,angles) 

Wrapped subroutines (so far)

  • makegriddh (expand spherical harmonic coefficients on a regular lat-lon grid)
  • shexpanddh (get spherical harmonics coefficients from a regular lat-lon grid)
  • shpowerspectrum (get power spectrum from spherical harmonics coefficients)
  • shcrosspowerspectrum (get cross power spectrum)
  • shadmitcorr (compute admittance and correlation)
  • shconfidence
  • shrotaterealcoef (rotate spherical harmonic coefficients by 3 angles)
For a detailed description of all the routines available in shtools see shtools.ipgp.fr .

Installation:

Prerequisites:

First you need to compile and install the shtools library. To do this, follow the instructions on its webpage (you'll also need fftw3 for this which is often available through the package manager of many distributions). I compiled with gfortran and changed the associated compiler flags in the shtools Makefile to:

ifeq ($(F95),gfortran)
# Default gfortran flags
F95FLAGS ?= -m64 -O3 -fPIC
MODFLAG = -Imodpath
endif 

 
The fPIC option is important and ensures compilation with position independent code.

Second you need to have scipy/numpy and especially the f2py compiler installed. You can test this by typing the command "f2py" in the console. It is often installed along numpy/scipy installations. Sometimes you have to set the correct path to it. You can find more info on how to use and install f2py on the google codes or scipy webpages.

Finally you can download the actual wrapper files here. wrapper.f90 contains the wrapper routines that need to be compiled with f2py. test.py is a small example script that shows how one can plot a map from spherical harmonic coefficients using shtools, numpy, matplotlib and basemap.

How to compile the wrapper

So far I only tested the compilation on my own linux machine. Once you have a working shtools library and also the f2py compiler running, try:


f2py -I"SHTOOLS_INCLUDE_DIR" -L"SHTOOLS_LIB_DIR" -lSHTOOLS2.8 -lfftw3 -lm --f90flags="-m64 -fPIC" --f77flags="-m64 -fPIC" -c -m shtools wrapper.f90 

this should create the file shtools.so that you can finally import into python.

The test script

the test script test.py shows an example on how to plot a map from spherical harmonics coefficients. You need shtools,numpy, matplotlib and basemap working. The script plots a function with sh-coefficients up to a certain degree. The result should look something like this:


7 comments:

  1. Hello Matthias,

    I begin to use SHTOOLS recently and I have a question concerning a function there. I sent emails to Dr. Wieczorek twice but couldn't receive an answer. May I ask you ? Sincerely

    ReplyDelete
  2. Hi Ali,

    yes of course. I am happy to help as far as I can!

    ReplyDelete
    Replies
    1. Hi Matthias,

      I try to compute Bouguer anomalies using SHTools. I try to follow Wieczorek (2007) and the SHTOOLS/examples/MarsCrustalThickness.
      First I compute the Free-air anomaly using eq. 22. This first step is successful. However, I find the Bouguer Correction output by CilmPlus quite weird.
      Could you help me in this issue? I can send the program that I wrote and input files.

      Sincerely,
      Ali

      Delete
  3. Hi Ali,

    I never used this routine so unfortunately I can't say much about it but what do you mean exactly with weird? Yes if you want send me your code (meschede at ipgp.fr). I only use the spherical harmonics transform and rotations, as well as the Wigner 3j functions and do all the rest in python. I have tested all these routines against 3rd party codes and they work very well! You should write to Mark immediately or to the SHTOOLS google group. I think they can help you more!

    ReplyDelete
    Replies
    1. Hi Matthias,

      By weird, I mean that there's a problem with scaling I presume. I wrote to Mark three times but he never answered back. It may be too easy a mistake I might be conducting but I am not sure where. Maybe I should try the google group.

      Thanks a lot.
      Ali

      Delete
  4. This comment has been removed by the author.

    ReplyDelete
    Replies
    1. Hi again,

      I posted the following message to the SHTOOLS email group. Unfortunately got no answer. I hope this paragrapg will explain the case better. Hope I can get a bit more help inthis regard.

      I began to study the topography - gravity relationships. In order to do this I will study Bouguer Anomaly (BA) admittance/coherence plots for the Eastern Mediterranean. I use the potential data from EGM2008 (Pavlis et al., 2008) and shape data from http://geodesy.curtin.edu.au/research/models/Earth2012/ : Earth2012.shape_bathy.SHCto2160.dat
      I am able to produce a free-air anomaly (FA) map with the following lines in the code (I multiply everything with gm/r^2 at the end):
      do l=0, lmaxp
      fa(:,l+1,1:l+1) = pot(:,l+1,1:l+1) * (l+1.0d0) * (R0_pot_Earth/r0)**l
      enddo
      ! Wieczorek, 2007, s. 9: zonal degree-2 term ... set to zero.
      fa(1:2,3,1:3)=0.0d0
      As the BA = FA - BC, I need to calculate the Bouguer Correction (including the terrain correction). I went through the code in examples/MarsCrustalThickness.f95, which briefly explains how to compute the BC using a function call to CilmPlus:
      call MakeGridDH(topo_grid, n_out, topo_c, lmax, norm = 1, &
      sampling = sampling, csphase = 1, lmax_calc = degmax)
      call CilmPlus(bc, topo_grid, degmax, nmax, mass, rref, rho_crust, gridtype, n = nlat)
      ba = (fa - bc) * gm / (R0_pot_Earth**2)
      However the resulting BA map is not correct. I suspect I make a mistake somewhere. Even though I remove the degree 1 topo coefficients from BC, the map does not become familiar.
      When I take a look in the MarsCrustalThickness.f95, I see that instead of calculating the FA, there potential coefficients are said to be down-propagated from R0_pot_Mars to R_Mars. That is different from FA calculation, for there is no multiplication of pot with (l+1) term, i.e. no differentiation. So I also tried that approach with the command -without understanding it-
      do l=1, degmax
      pot(1:2,l+1,1:l+1) = pot(1:2,l+1,1:l+1) * (R0_pot_Earth/r0)*l ! 1.0d5
      enddo
      But neither did this method worked for me. I should say that I can reproduce the Martian anomalies perfectly well (with comparing figures from Wieczorek, 2007).
      I would appreciate help in all the matter aforementioned. If requested I can post the entire program here.
      With kindest regards,
      Ali Ozbakir.

      Delete