From 2e602a97bba182c350bcdb51d9eca539a2adf156 Mon Sep 17 00:00:00 2001 From: Martin Jucker Date: Wed, 9 Sep 2015 12:23:03 -0400 Subject: [PATCH] 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 --- basic.py | 37 ++++++++++++++++++++++++++++--------- examples/example_flat.py | 3 ++- examples/example_ocean.py | 5 +++-- examples/example_pole.py | 6 +++--- examples/example_sphere.py | 3 ++- grids.py | 25 ++++++++++++++----------- testing/test_functions.py | 29 ++++++++++++++++------------- 7 files changed, 68 insertions(+), 40 deletions(-) diff --git a/basic.py b/basic.py index f1a08ac..b98aabd 100644 --- a/basic.py +++ b/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 : diff --git a/examples/example_flat.py b/examples/example_flat.py index 708270c..3686af4 100644 --- a/examples/example_flat.py +++ b/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 * diff --git a/examples/example_ocean.py b/examples/example_ocean.py index c8e1613..acdca28 100644 --- a/examples/example_ocean.py +++ b/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() diff --git a/examples/example_pole.py b/examples/example_pole.py index 75d91f9..61b9c46 100644 --- a/examples/example_pole.py +++ b/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() diff --git a/examples/example_sphere.py b/examples/example_sphere.py index 4a91b2b..bb8cf93 100644 --- a/examples/example_sphere.py +++ b/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 * diff --git a/grids.py b/grids.py index 364b584..d962be6 100644 --- a/grids.py +++ b/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) diff --git a/testing/test_functions.py b/testing/test_functions.py index 5845603..50a7943 100644 --- a/testing/test_functions.py +++ b/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'