More general coordinate handling

- logCoords and basis are now arrays to allow for multiple logarithmic
coordinates instead of only one.
- there are now additional inputs to LoadData (and thus also
TransformCoords): reverseCoords and revCenter. These are arrays that
tell pv_atmos which coordinates should be reversed (such as e.g.
pressure), and around which value

- The example files as well as the test_functions file now ask for the
path to pv_atmos instead of assuming something.
- Example and test files are modified to take into account the changes
to pv_atmos mentioned above

Author: mjucker

basic.py

@@ -30,7 +30,7 @@ def ConvertLogCoordString(pString, basis=1e3):
     return expression
 
 # do the coordinate conversion inside a Calculator
-def Cart2Log(src=GetActiveSource(), ratios=[1,1,1], logCoords=[2], basis=[1e3]):
+def Cart2Log(src=GetActiveSource(), ratios=[1,1,1], logCoords=[], basis=[]):
     """Convert between logarithmic and linear coordinates. Also applies aspect ratio correction.
 
     Adds a Calculator filter to the pipeline
@@ -86,16 +86,33 @@ def GridAspectRatio(ratios, src=GetActiveSource()):
     ratios -- 2- or 3-vector with multiplicative factors for each spatial coordinate
     """
     calc=Calculator(src)
-    try:
-        calc.Function = 'iHat*'+str(ratios[0])+'*coordsX + jHat*'+str(ratios[1])+'*coordsY + kHat*'+str(ratios[2])+'*coordsZ'
-    except:
-        calc.Function = 'iHat*'+str(ratios[0])+'*coordsX + jHat*'+str(ratios[1])+'*coordsY'
+    calc.Function = 'iHat*'+str(ratios[0])+'*coordsX + jHat*'+str(ratios[1])+'*coordsY + kHat*'+str(ratios[2])+'*coordsZ'
     calc.CoordinateResults = 1
     return calc
 
 # transform coordinates: logarithmic, aspect ratio
-def TransformCoords(src=GetActiveSource(), aspectRatios=[1,1,1], logCoords=[2], basis=[1e3]):
+def TransformCoords(src=GetActiveSource(), aspectRatios=[1,1,1], logCoords=[], basis=[], reverseCoords=[], revCenter=[]):
     """Transform the coordinates depending on whether or not there are logarithmic coordinates"""
+    if len(reverseCoords)>0:
+        nVec = ['iHat*','jHat*','kHat*']
+        nCoor= ['X','Y','Z']
+        revCoor = Calculator(src)
+        rFun = ''
+        for dim in range(3):
+            if dim in reverseCoords or -dim in reverseCoords:
+                for d in range(len(reverseCoords)):
+                    if dim == abs(reverseCoords[d]):
+                        rd = d
+                if reverseCoords[rd]<0:
+                    coorSign = '+'
+                else:
+                    coorSign = '-'
+                rFun += ' +'+nVec[dim]+'('+str(revCenter[rd])+coorSign+'coords'+nCoor[dim]+')'
+            else:
+                rFun += ' +'+nVec[dim]+'coords'+nCoor[dim]
+        revCoor.Function = rFun[2:]
+        revCoor.CoordinateResults = 1
+        src = revCoor
     if len(logCoords)>0 :
         transCoor = Cart2Log(src=src,ratios=aspectRatios,logCoords=logCoords,basis=basis)
     else:
@@ -111,7 +128,7 @@ def MakeSelectable(src=GetActiveSource()):
 
 ######### read in data, redefine pressure coordinates and change aspect ratio ###############
 
-def LoadData( fileName, ncDims=['lon','lat','pfull'], aspectRatios=[1,1,1], logCoords=[2], basis=[1e3], replaceNaN=True ):
+def LoadData( fileName, ncDims=['lon','lat','pfull'], aspectRatios=[1,1,1], logCoords=[], basis=[], reverseCoords=[], revCenter=[], replaceNaN=True ):
     """Load netCDF file, convert coordinates into useful aspect ratio.
 
     Adds file output_nc, and Calculator LogP or Calculator AspRat to the pipeline
@@ -120,8 +137,10 @@ def LoadData( fileName, ncDims=['lon','lat','pfull'], aspectRatios=[1,1,1], logC
         fileName      -- full path and file name of data to be read
         ncDims        -- names of the dimensions within the netCDF file. Time should be excluded. Ordering [x,y,z]
         aspectRatios  -- how to scale coordinates [xscale,yscale,zscale]. Z coordinate is scaled after applying log10 for logarithmic axes
-        logCoords     -- index/indices of dimension(s) to be logarithmic
+        logCoords     -- index/indices of dimension(s) to be logarithmic, set to [] if no log coordinates
         basis         -- basis to normalize argument to logarithm (ie defines origin). List of same length as logCoords
+        reverseCoords -- index/indices of dimension(s) to be reversed, set to [] if none to be reversed
+        revCenter     -- center of reversal if reverseCoords is not empty. List of same length as logCoords
         replaceNaN    -- whether or not to replace the FillValue with NaNs
     OUTPUTS:
         output_nc     -- netCDF reader object with the file data as read
@@ -144,7 +163,7 @@ def LoadData( fileName, ncDims=['lon','lat','pfull'], aspectRatios=[1,1,1], logC
     MakeSelectable()
     RenameSource(fileName,output_nc)
 
-    transCoor = TransformCoords(src=output_nc,aspectRatios=aspectRatios,logCoords=logCoords,basis=basis)
+    transCoor = TransformCoords(src=output_nc,aspectRatios=aspectRatios,logCoords=logCoords,basis=basis,reverseCoords=reverseCoords,revCenter=revCenter)
     MakeSelectable()
 
     if len(logCoords)>0 :

examples/example_flat.py

@@ -12,7 +12,8 @@
 try:
     pvAtmosPath
 except:
-    pvAtmosPath = '../'
+    pvAtmosPath = raw_input("Please provide the path where pv_atmos lives: ")
+    pvAtmosPath = pvAtmosPath+'/'
 try: #is pv_atmos installed?
     from pv_atmos.basic import *
     from pv_atmos.grids import *

examples/example_ocean.py

@@ -12,7 +12,8 @@
 try:
     pvAtmosPath
 except:
-    pvAtmosPath = '../'
+    pvAtmosPath = raw_input("Please provide the path where pv_atmos lives: ")
+    pvAtmosPath = pvAtmosPath+'/'
 try: #is pv_atmos installed?
     from pv_atmos.basic import *
     from pv_atmos.grids import *
@@ -52,7 +53,7 @@
 
 
 # bathymetry
-(depth_out,depth_coor)=LoadData(oceanPath+topoFile,ncDims=topoDims,logCoords=logCoord )
+(depth_out,depth_coor)=LoadData(oceanPath+topoFile,ncDims=topoDims,logCoords=logCoord,replaceNaN=False )
 # get the bounds of the topography. This is important if bathymetry and data files have different origins
 topoBds = depth_out.GetDataInformation().GetBounds()
 

examples/example_pole.py

@@ -12,7 +12,8 @@
 try:
     pvAtmosPath
 except:
-    pvAtmosPath = '../'
+    pvAtmosPath = raw_input("Please provide the path where pv_atmos lives: ")
+    pvAtmosPath = pvAtmosPath+'/'
 try: #is pv_atmos installed?
     from pv_atmos.basic import *
     from pv_atmos.grids import *
@@ -45,7 +46,6 @@
 
 # then, the file containing geopotential height data
 dataFile = 'ECMWF_19790223.nc'
-dataFile = 'tmp.nc'
 dataDims = ['longitude','latitude','level']
 # the values we will be interested in
 dataName = 'height'
@@ -68,7 +68,7 @@
 ## Earth
 
 # get bathymetry: this is exactly the same as in example_ocean.py
-(depth_out,depth_coor)=LoadData(topoFileName,ncDims=topoDims,logCoords=topoLogCoord )
+(depth_out,depth_coor)=LoadData(topoFileName,ncDims=topoDims,logCoords=topoLogCoord, replaceNaN = False )
 # the bathymetry file is in cell data, need to convert to point data
 c2p=CellDatatoPointData(depth_coor)
 MakeSelectable()

examples/example_sphere.py

@@ -12,7 +12,8 @@
 try:
     pvAtmosPath
 except:
-    pvAtmosPath = '../'
+    pvAtmosPath = raw_input("Please provide the path where pv_atmos lives: ")
+    pvAtmosPath = pvAtmosPath+'/'
 try: #is pv_atmos installed?
     from pv_atmos.basic import *
     from pv_atmos.grids import *

grids.py

@@ -396,22 +396,25 @@ def SphericalShells(radius=1, ratios=[1,1,1], logCoord=[2], basis=[1e3], src=Get
     return Planes
 
 ###### add full set of grids and lables in rectangular geometry ############################
-def AddGrid(xlevels=[0,90,180,270], ylevels=[-60,-30,0,30,60], zlevels=[100,10,1,0.1], bounds=[0.0,360.0,-90.0,90.0,1e3,0.1], ratios=[1,1,1], logCoord=[2], basis=[1e3], AxisNames=["lon","lat","pressure [hPa]"], AxisColor=[0,0,0], AxisWidth=2.0,LabelSize=5.0):
+def AddGrid(xlevels=[0,90,180,270], ylevels=[-60,-30,0,30,60], zlevels=[100,10,1,0.1], bounds=[0.0,360.0,-90.0,90.0,1e3,0.1], ratios=[1,1,1], logCoord=[], basis=[], reverseCoords=[], revCenter=[], AxisNames=["lon","lat","pressure [hPa]"], AxisColor=[0,0,0], AxisWidth=2.0,LabelSize=5.0):
     """Add a full grid with grid lines and labels.
         
         Adds as many X,Y,Z grid lines as needed. This function adds a lot of objects and filters to the pipeline, and should probably only be used once the visualization itself is finished. This function can be called even if there is no data loaded.
 
         INPUTS:
-            xlevels   -- vector with X grid positions
-            ylevels   -- vector with Y grid positions
-            zlevels   -- vector with Z grid levels
-            bounds    -- grid bounds
-            ratios    -- axes ratios
-            basis     -- basis (surface) pressure
-            AxisNames -- names of x,y,z axis
-            AxisColor -- color of lines in RGB
-            AxisWidth -- width of lines
-            LabelSize -- Size of the label text
+            xlevels       -- vector with X grid positions
+            ylevels       -- vector with Y grid positions
+            zlevels       -- vector with Z grid levels
+            bounds        -- grid bounds
+            ratios        -- axes ratios
+            logCoord      -- coordinate index/indices for logarithmic axes
+            basis         -- basis for logarithmic coordinates logCoords
+            reverseCoords -- coordinate index/indices for axes to be reversed
+            revCenter     -- center around which to reverse reverseCoords
+            AxisNames     -- names of x,y,z axis
+            AxisColor     -- color of lines in RGB
+            AxisWidth     -- width of lines
+            LabelSize     -- Size of the label text
         """
     (Xmin,Xmax,Ymin,Ymax,Zmin,Zmax) = BoundAspectRatio(bounds, ratios, logCoord, basis)
     absLsz = abs(LabelSize)

testing/test_functions.py

@@ -4,22 +4,25 @@
 import numpy as np
 import os
 
-#adjust the path to pv-atmos here:
-try:
-	pvAtmosPath=os.path.dirname(__file__) + '/../'
-except: #if run in interactive mode
-	pvAtmosPath='../'
 try: from atmos_basic import *
-except: execfile(pvAtmosPath + 'atmos_basic.py')
+except:
+    pvAtmosPath = raw_input("Please provide the path where pv_atmos lives: ")
+    pvAtmosPath = pvAtmosPath+'/'
+    print 'Using '+pvAtmosPath
+    execfile(pvAtmosPath + 'basic.py')
 try: from atmos_grids import *
-except: execfile(pvAtmosPath + 'atmos_grids.py')
+except: execfile(pvAtmosPath + 'grids.py')
 
 
 class TestSequenceFunctions(unittest.TestCase):
     def setUp(self):
         self.p = random.uniform(1e-5,1e3) # choose pressure surface randomly
         self.rat = [random.uniform(0,1),random.uniform(0,1),random.uniform(0,1)] # choose box aspect ration randomly
-        self.pCoord = round(random.uniform(0,1)) # choose whether or not data uses pressure coordinates
+        # choose whether or not data uses pressure coordinates
+        if round(random.uniform(0,1)) == 1:
+            self.pCoord = [int(round(random.uniform(0,2)))] # which dimension is logarithmic?
+        else:
+            self.pCoord = []
         self.np = round(random.uniform(1,5)) # number of pressure levels for grid
         self.nx = round(random.uniform(1,5)) # number of lon gridlines for grid
         self.ny = round(random.uniform(1,5)) # number of lat gridlines for grid
@@ -35,21 +38,21 @@ def setUp(self):
     def test_basic(self):
         # load data
         fileName = pvAtmosPath+'examples/uv_daily.nc'
-        (output_nc,CorrZ,Coor,AspRat) = loadData(fileName, ['pfull','lat','lon'], self.pCoord, self.rat, math.ceil(self.p))
+        (output_nc,AspRat) = LoadData(fileName, ncDims=['lon','lat','pfull'], aspectRatios=self.rat, logCoords=self.pCoord, basis=[math.ceil(self.p)])
         self.assertEqual(output_nc.OutputType,'Unstructured')
-        (W,norm,clipS,clipN)=CartWind2Atmos(src=AspRat,ratios=self.rat)
+        (W,norm,clipS,clipN)=CartWind2Sphere(src=AspRat,ratios=self.rat)
         self.assertEqual(clipS.ClipType.Origin[1],-80.0*self.rat[1])
         self.assertEqual(clipN.ClipType.Origin[1], 80.0*self.rat[1])
 
     def test_allGrids(self):
-        AddGrid(press=np.arange(self.p,self.top,-self.p/self.np), lats=np.arange(self.mx[0],self.mx[1],(self.mx[1]-self.mx[0])/self.nx), lons=np.arange(self.my[0],self.my[1],(self.my[1]-self.my[0])/self.ny), bounds=[self.my[0],self.my[1],self.mx[0],self.mx[1],self.p,self.top], ratios=self.rat, basis=math.ceil(self.p), AxisColor=[0,0,0], AxisWidth=2.0,LabelSize=5.0)
+        AddGrid(zlevels=np.arange(self.p,self.top,-self.p/self.np), ylevels=np.arange(self.mx[0],self.mx[1],(self.mx[1]-self.mx[0])/self.nx), xlevels=np.arange(self.my[0],self.my[1],(self.my[1]-self.my[0])/self.ny), bounds=[self.my[0],self.my[1],self.mx[0],self.mx[1],self.p,self.top], ratios=self.rat, basis=[math.ceil(self.p)], AxisColor=[0,0,0], AxisWidth=2.0,LabelSize=5.0)
         print 'Check pipeline for grid'
 
     def test_shells(self):
         # load data
         fileName = pvAtmosPath+'examples/uv_daily.nc'
-        (output_nc,CorrZ,Coor,AspRat) = loadData(fileName, ['pfull','lat','lon'], 1, self.rat, math.ceil(self.p))
-        AtmosShells(radius=self.rad,ratios=self.rat,basis=math.ceil(self.p),src=AspRat,shellValues=np.arange(self.p,self.top,-self.p/self.np),labels=self.pCoord,labelPosition=[self.my[0],self.mx[0]],waterMark='WaterMark',markPosition=[self.my[1],self.mx[1]])
+        (output_nc,AspRat) = LoadData(fileName, ['pfull','lat','lon'], self.rat, [1], [math.ceil(self.p)])
+        SphericalShells(radius=self.rad,ratios=self.rat,basis=[math.ceil(self.p)],src=AspRat,shellValues=np.arange(self.p,self.top,-self.p/self.np),labels=self.pCoord,labelPosition=[self.my[0],self.mx[0]],waterMark='WaterMark',markPosition=[self.my[1],self.mx[1]])
         print 'Check pipeline for spherical shells'