These examples show how to use some of the available components of the python module. They uses the same dataset as for Section 3.3, “A Fortran 90 example”. If you configurations allows it, you can run theses examples typing for example this from the installation package directory:
<user> cd example && make python1
<user> make install
This first example is similar to Section 3.3, “A Fortran 90 example”: a PCA/MSSA is performed on the whole dataset, and the first oscillation (pair of MSSA modes) is reconstructed using phases composites (16 phases).
Example 2. Python example
# File: example1.py # # This file is part of the SpanLib library. # Copyright (C) 2006 Charles Doutiraux, Stephane Raynaud # Contact: stephane dot raynaud at gmail dot com # # This library is free software; you can redistribute it and/or # modify it under the terms of the GNU Lesser General Public # License as published by the Free Software Foundation; either # version 2.1 of the License, or (at your option) any later version. # # This library is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public # License along with this library; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA ################################################################### # In this example, we perform a pre-PCA to reduce the number of # d-o-f, then we perform an MSSA to extract the first oscillation # (first par of modes). Finally, we compute phase composites # to represent the oscillation over its cycle. ################################################################### print "#################################################" print "# PCA+MSSA+Phase_composites of 1st oscillation. #" print "# Then reconstructions, and plots: #" print "# - 1 phase over 2 #" print "# - 1 time series #" print "#################################################" # Needed modules import sys import cdms import spanlib import vcs import MV import Numeric import cdutil import genutil # We tell cdms that we have longitude, latitude and time cdms.axis.latitude_aliases.append('Y') cdms.axis.longitude_aliases.append('X') cdms.axis.time_aliases.append('T') # Simply open the netcdf file print "Open file" f=cdms.open('../example/data2.cdf') # Retrieve data print "Read the whole dataset" s=f('ssta',time=slice(0,120)) # Create the analysis object print "Creating SpAn object" SP=spanlib.SpAn(s) # Perform a preliminary PCA+MSSA # (equivalent to simple use SP.mssa(pca=True) later) print "PCA..." eof,pc,ev = SP.pca() # MSSA on PCA results print 'MSSA...' steof,stpc,stev = SP.mssa() # Phase composites of first two MSSA modes print 'Phase composites...' out = SP.reconstruct(phases=True,nphases=16,end=2) # Plot 1 phase over two, then a time series print "Now, plot!" x=vcs.init() for i in range(0,out.shape[0],2): x.plot(out[i],title="Phase composites of the first MSSA oscillation") raw_input('map out %i/%i ok?' % ( i+1 , out.shape[0])) x.clear() x.plot(out[:,30,80],title="Cycle of the ocillation") raw_input('Time series at center of bassin ok?') x.clear()
In this second example, two different areas are analysed at the same time. This is intented to mimics the used of two different datasets that are stacked, before being analysed. Such approach (see for example Raynaud et al (2006)) allows to find modes of variability in a arbitrary number of variables, provided you are careful with units (performing appropriate normalisations).
Example 3. Python example
# File: example2.py # # This file is part of the SpanLib library. # Copyright (C) 2006 Charles Doutiraux, Stephane Raynaud # Contact: stephane dot raynaud at gmail dot com # # This library is free software; you can redistribute it and/or # modify it under the terms of the GNU Lesser General Public # License as published by the Free Software Foundation; either # version 2.1 of the License, or (at your option) any later version. # # This library is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public # License along with this library; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA ################################################################### # In this example, we analyse two different areas at the same time. # You can do the same with two completely different datasets, # except that they must have the same temporal grid. ################################################################### print "##############################################" print "# PCA+MSSA applied on two different regions. #" print "# Then reconstructions and plots. #" print "##############################################" # Needed modules import cdms import spanlib import MV import vcs # We tell cdms that we have longitude, latitude and time cdms.axis.latitude_aliases.append('Y') cdms.axis.longitude_aliases.append('X') cdms.axis.time_aliases.append('T') # Simply open the netcdf file print "Open file" f=cdms.open('data2.cdf') # Get our two datasets print "Read two different regions" s2=f('ssta',latitude=(-10,10),longitude=(110,180)) s1=f('ssta',latitude=(-15,15),longitude=(210,250)) # Stack the two dataset to have only one dataset print "Stacking data" res = spanlib.stackData(s1,s2) # Create the analysis object print "Creating SpAn object" SP=spanlib.SpAn(MV.array(res[0]),weights=MV.array(res[1])) # Perform a preliminary PCA # (optional step since done by default with mssa) print "PCA+MSSA..." steof,stpc,stev = SP.mssa(pca=True) # Recontructed the filtered field print "Reconstructing selected modes" ffrec = SP.reconstruct() # Unstacking print "Unstaking data" out = spanlib.unStackData(ffrec,res[1],res[2],res[3]) # Plot a timeseries taken from our two # recontructed datasets print "Time series for the two filtered regions" x=vcs.init() x.plot(out[1][:,5,5]) x.plot(out[0][:,5,5]) raw_input('ok?')