From e26251f2204be094456c1db78d2edc9224c78ba6 Mon Sep 17 00:00:00 2001 From: Soren Rasmussen Date: Tue, 23 Apr 2024 12:41:11 -0600 Subject: [PATCH] Update to README on how to setup PYTHONPATH and what scripts use deprecated units package. Small fixes to other python scripts so they run. README updated to say that the Nio package is not installed as a dependency here --- helpers/README.md | 22 ++++++++++-- helpers/ccsm/config.py | 20 +++++------ helpers/cmip/io_routines.py | 30 ++++++++-------- helpers/erai/config.py | 5 ++- helpers/gen_bc.py | 24 ++++++------- helpers/gen_ideal.py | 45 ++++++++++++----------- helpers/gen_init.py | 26 +++++++------- helpers/gen_init_ideal.py | 8 ++--- helpers/wrf/ideal2icar.py | 32 ++++++++--------- helpers/wrf/wind_compare.py | 71 ++++++++++++++++++------------------- helpers/wrf/wrf_vars.py | 10 +++--- 11 files changed, 151 insertions(+), 142 deletions(-) diff --git a/helpers/README.md b/helpers/README.md index 238c73c1..5aecafcc 100644 --- a/helpers/README.md +++ b/helpers/README.md @@ -1,15 +1,31 @@ # Install Python Dependencies The following instructions and dependecy files work for the core ICAR scripts. -Tools in `make_domain.py` and ccsm, cesm, cmip, erai, and wrf directories will require Bunch and mygis packages as well. +Tools in `make_domain.py` and ccsm, cesm, cmip, erai, and wrf directories will require the mygis packages as well. The Python script `ideal_linear.py` will require Nio to be installed with `pip install nio`. -## Install With Conda +## Setup Environment +### Install With Conda ```bash $ conda env create -f environment.yml --prefix /path/to/install/icar_env $ conda activate icar_env + +set PYTHONPATH, this will be saved for future use by the environment +$ conda env config vars set PYTHONPATH=$(pwd)/lib:$PYTHONPATH + +reactivate environment +$ conda activate icar_env ``` -## Install With Pip +### Install With Pip ```bash $ pip install -r requirements.txt ``` +Make sure the `lib` directory is in the `PYTHONPATH`, add to `.bashrc` or other startup files for repeat use. +```bash +$ export PYTHONPATH=$(pwd)/lib:$PYTHONPATH +``` + + +## Deprecated Scripts +The units package used in the `cmip/convert.py, cesm/bias_correct.py` and `gen_sounding.py` scripts has been deprecated. +The [Nio package](https://www.pyngl.ucar.edu/Nio.shtml) used in `create_geo_testfiles.py` is not installed as a dependency. diff --git a/helpers/ccsm/config.py b/helpers/ccsm/config.py index b423db2b..d48d6fef 100644 --- a/helpers/ccsm/config.py +++ b/helpers/ccsm/config.py @@ -15,28 +15,28 @@ def set_bounds(info): ccsm_file=atm_file.replace("_Y_","2006").replace("_M_","01").replace("_D_","01").replace("_VAR_","hus") ccsm_file=glob.glob(ccsm_file)[0] varlist=["lat","lon"] - + lat=io.read_nc(ccsm_file,varlist[0]).data lon=io.read_nc(ccsm_file,varlist[1]).data-360 - + info.xmin=np.where(lon>=info.lon[0])[0][0] info.xmax=np.where(lon<=info.lon[1])[0][-1]+1 info.ymin=np.where(lat>=info.lat[0])[0][0] info.ymax=np.where(lat<=info.lat[1])[0][-1]+1 - + lon,lat=np.meshgrid(lon[info.xmin:info.xmax],lat[info.ymin:info.ymax]) info.lat_data=lat info.lon_data=lon - + def make_timelist(info,hrs=6.0): dt=datetime.timedelta(hrs/24) - info.ntimes=np.int(np.round((info.end_date-info.start_date).total_seconds()/60./60./hrs)) + info.ntimes=np.int64(np.round((info.end_date-info.start_date).total_seconds()/60./60./hrs)) info.times=[info.start_date+dt*i for i in range(info.ntimes)] def update_info(info): make_timelist(info) set_bounds(info) - + def parse(): parser= argparse.ArgumentParser(description='Convert CCSM files to ICAR input forcing files') @@ -51,24 +51,24 @@ def parse(): parser.add_argument('sfcdir', nargs="?",action='store',help="CCSM surface data file location", default="ccsm_sfc/") parser.add_argument('atmfile', nargs="?",action='store',help="CCSM atmospheric files", default="_VAR__6hrLev_CCSM4_rcp85_r6i1p1__Y__M__D_00-*.nc") parser.add_argument('sfcfile', nargs="?",action='store',help="CCSM surface files", default="_VAR__3hr_CCSM4_rcp85_r6i1p1__Y__M__D_0000-*.nc") - + parser.add_argument('-v', '--version',action='version', version='CCSM2ICAR v'+version) parser.add_argument ('--verbose', action='store_true', default=False, help='verbose output', dest='verbose') args = parser.parse_args() - + date0=args.start_date.split("-") start_date=datetime.datetime(int(date0[0]),int(date0[1]),int(date0[2])) date0=args.end_date.split("-") end_date=datetime.datetime(int(date0[0]),int(date0[1]),int(date0[2])) - + info=Bunch(lat=[float(args.lat_s),float(args.lat_n)], lon=[float(args.lon_w),float(args.lon_e)], start_date=start_date, end_date=end_date, atmdir=args.dir+args.atmdir, sfcdir=args.dir+args.sfcdir, atmfile=args.atmfile, sfcfile=args.sfcfile, version=version) - + return info diff --git a/helpers/cmip/io_routines.py b/helpers/cmip/io_routines.py index 8d5d9622..a6aa0227 100644 --- a/helpers/cmip/io_routines.py +++ b/helpers/cmip/io_routines.py @@ -17,9 +17,9 @@ def read_nc(filename,var="data",proj=None,returnNCvar=False): data:raw data as an array proj:string representation of the projection information atts:data attribute dictionary (if any) - if (returnNCvar==True) then the netCDF4 file is note closed and the netCDF4 - representation of the variable is returned instead of being read into - memory immediately. + if (returnNCvar==True) then the netCDF4 file is note closed and the netCDF4 + representation of the variable is returned instead of being read into + memory immediately. ''' d=netCDF4.Dataset(filename, mode='r',format="nc") outputdata=None @@ -39,8 +39,8 @@ def read_nc(filename,var="data",proj=None,returnNCvar=False): if proj!=None: projection=d.variables[proj] outputproj=str(projection) - - + + if returnNCvar: return Bunch(data=outputdata,proj=outputproj,ncfile=d,atts=attributes) d.close() @@ -53,7 +53,7 @@ def find_atm_file(time,varname,info): file_base= file_base.replace("_Y_",str(time.year)) file_base= file_base.replace("_EXP_",info.experiment) atm_file = file_base.replace("_ENS_",info.ensemble) - + print(atm_file) filelist = glob.glob(atm_file) filelist.sort() @@ -67,7 +67,7 @@ def find_sst_file(time,info): file_base= file_base.replace("_Y_",str(time.year)) file_base= file_base.replace("_EXP_",info.experiment) sst_file = file_base.replace("_ENS_",info.ensemble) - + print(sst_file) filelist=glob.glob(sst_file) filelist.sort() @@ -76,7 +76,7 @@ def find_sst_file(time,info): def load_atm(time,info): """Load atmospheric variable from a netcdf file""" - + outputdata=Bunch() for s,v in zip(icar_atm_var,atmvarlist): @@ -98,24 +98,24 @@ def load_atm(time,info): outputdata[varname]=np.concatenate([outputdata[varname],newdata]) else: outputdata[varname]=newdata - + varname="p" newdata=info.read_pressure(atmfile)[:,:,info.ymin:info.ymax,info.xmin:info.xmax] if varname in outputdata: outputdata[varname]=np.concatenate([outputdata[varname],newdata]) else: outputdata[varname]=newdata - + outputdata.ntimes = outputdata.p.shape[0] - + # outputdata.times=info.read_time(atmfile) try: calendar = mygis.read_attr(atmfile_list[0], "calendar", varname="time") - except KeyError,IndexError: + except (KeyError,IndexError): calendar = None - + outputdata.calendar = calendar - + return outputdata def load_sfc(time,info): @@ -137,5 +137,3 @@ def load_data(time,info): atm=load_atm(time,info) sfc=load_sfc(time,info) return Bunch(sfc=sfc,atm=atm) - - diff --git a/helpers/erai/config.py b/helpers/erai/config.py index 3eb8a96a..535cf2b3 100644 --- a/helpers/erai/config.py +++ b/helpers/erai/config.py @@ -3,9 +3,8 @@ import numpy as np -from bunch import Bunch import sys -sys.path.append('../lib') # required to find mygis.py. Update as required. +from bunch import Bunch import mygis import io_routines as io @@ -53,7 +52,7 @@ def set_bounds(info): def make_timelist(info): hrs=6.0 dt=datetime.timedelta(hrs/24) - info.ntimes=np.int(np.round((info.end_date-info.start_date).total_seconds()/60./60./hrs)) + info.ntimes=np.int64(np.round((info.end_date-info.start_date).total_seconds()/60./60./hrs)) info.times=[info.start_date+dt*i for i in range(info.ntimes)] def update_info(info): diff --git a/helpers/gen_bc.py b/helpers/gen_bc.py index a07070be..7ce87c52 100755 --- a/helpers/gen_bc.py +++ b/helpers/gen_bc.py @@ -26,8 +26,8 @@ def update_base(base,filename,nz): def main(): filename="bc" nx,ny,nz,nt=(20.,20,10,24) - dims=[nt,nz,ny,nx] - + dims=[nt,nz,ny,int(nx)] + lonmin=-110.0; lonmax=-100.0; dlon=(lonmax-lonmin)/nx latmin=35.0; latmax=45.0; dlat=(latmax-latmin)/ny @@ -41,32 +41,32 @@ def main(): update_base(base,"sounding.txt",nz) nz=base.th.size dims=[nt,nz,ny,nx] - + u=np.zeros(dims,dtype="f")+base.u w=np.zeros(dims,dtype="f")+base.w v=np.zeros(dims,dtype="f")+base.v qv=np.zeros(dims,dtype="f")+base.qv qc=np.zeros(dims,dtype="f")+base.qc coscurve=np.cos(np.arange(dims[2])/dims[2]*2*np.pi+np.pi)+1 - hgt=(coscurve*1000).reshape((1,nx)).repeat(ny,axis=0) - + hgt=(coscurve*1000).reshape((1,int(nx))).repeat(ny,axis=0) + lon=np.arange(lonmin,lonmax,dlon) lat=np.arange(latmin,latmax,dlat) lon,lat=np.meshgrid(lon,lat) - + dz=np.zeros(dims)+base.dz - z=np.zeros(dims,dtype="f")+base.z.reshape((1,nz,1,1))+hgt.reshape((1,1,ny,nx)) - + z=np.zeros(dims,dtype="f")+base.z.reshape((1,nz,1,1))+hgt.reshape((1,1,ny,int(nx))) + layer1=(dz[0,0,:,:]/2) z[0,0,:,:]+=layer1 for i in range(1,int(nz)): z[:,i,:,:]=z[:,i-1,:,:]+(dz[:,i-1,:,:]+dz[:,i,:,:])/2.0 - + p=np.zeros(dims,dtype="f")+base.p adjust_p(p,0.0,z) th=np.zeros(dims,dtype="f")+base.th - - + + d4dname=("t","z","y","x") d3dname=("z","y","x") d2dname=("y","x") @@ -86,7 +86,7 @@ def main(): if fileexists: print("Removing : "+fileexists[0]) os.remove(fileexists[0]) - + io.write(filename, u,varname="U", dims=d4dname,dtype="f",attributes=dict(units="m/s",description="Horizontal (x) wind speed"), extravars=othervars) diff --git a/helpers/gen_ideal.py b/helpers/gen_ideal.py index 5ebc3a58..2369b45a 100755 --- a/helpers/gen_ideal.py +++ b/helpers/gen_ideal.py @@ -30,7 +30,7 @@ def adjust_p(p,h,dz): def update_base(base,filename,nz): """update the base information using data from a sounding file - + filename should be a space delimited text file with 3 columns height [m], potential temperature [K], and specific humidity [g/kg]""" print("Using Sounding from : "+filename) @@ -43,7 +43,7 @@ def update_base(base,filename,nz): def build_topography(experiment,dims): """create the topography to be used for a given case study""" - + # experiment D1 D2 D3 # h height of the hill 1800 1400 1040 [meters] # sigma half-width 60 40 3.1 [grid cells] @@ -60,8 +60,8 @@ def build_topography(experiment,dims): x = np.linspace(0,Lx,Nx) # % distance array (m) zo = [1700.0,2000.0,2200.0][experiment] # mountain base height (m) NOT REALLY USED CORRECTLY YET, NEED to truncate the sounding... zo = 0.0 - - + + zs=hm/(1.0+((x/dx-xm)/am)**2.) # zs-=zs[0] zs=zs.reshape((1,dims[3])).repeat(dims[2],axis=0) @@ -71,8 +71,8 @@ def main(): filename="ideal_{}_{}".format(case_study,int(wind_speed)) print(filename) nx,nz,ny=master_dims[case_study] - dims=[1,nz,ny,nx] - + dims=[1,int(nz),int(ny),int(nx)] + # this is just arbitrary for now dlon=dx/111.1 dlat=dx/111.1 @@ -89,7 +89,7 @@ def main(): update_base(base,"sounding.txt",nz) nz=base.th.size dims=[1,nz,ny,nx] - + udims=copy(dims) udims[-1]+=1 vdims=copy(dims) @@ -99,40 +99,40 @@ def main(): v=np.zeros(vdims,dtype="f")+base.v qv=np.zeros(dims,dtype="f")+base.qv qc=np.zeros(dims,dtype="f")+base.qc - + # simple topography = a cosine # coscurve=np.cos(np.arange(dims[3])/dims[3]*2*np.pi+np.pi)+1 # hgt=(coscurve*1000).reshape((1,nx)).repeat(ny,axis=0) hgt=build_topography(case_study,dims) - - lon=np.arange(lonmin,lonmax,dlon)[:nx] + + lon=np.arange(lonmin,lonmax,dlon)[:int(nx)] lat=np.arange(latmin,latmax,dlat)[:ny] lon,lat=np.meshgrid(lon,lat) - ulon=np.arange(lonmin-dlon/2,lonmax+dlon/2,dlon)[:nx+1] + ulon=np.arange(lonmin-dlon/2,lonmax+dlon/2,dlon)[:int(nx)+1] ulat=np.arange(latmin,latmax,dlat)[:ny] ulon,ulat=np.meshgrid(ulon,ulat) - vlon=np.arange(lonmin,lonmax,dlon)[:nx] + vlon=np.arange(lonmin,lonmax,dlon)[:int(nx)] vlat=np.arange(latmin-dlat/2,latmax+dlat/2,dlat)[:ny+1] vlon,vlat=np.meshgrid(vlon,vlat) - + dz=np.zeros(dims)+base.dz - z=np.zeros(dims,dtype="f")+base.z.reshape((1,nz,1,1))+hgt.reshape((1,1,ny,nx)) - + z=np.zeros(dims,dtype="f")+base.z.reshape((1,nz,1,1))+hgt.reshape((1,1,ny,int(nx))) + layer1=(dz[0,:,:]/2) z[0,:,:]+=layer1 for i in range(1,int(nz)): z[:,i,:,:]=z[:,i-1,:,:]+(dz[:,i-1,:,:]+dz[:,i,:,:])/2.0 - + p=np.zeros(dims,dtype="f")+base.p adjust_p(p,0.0,z) th=np.zeros(dims,dtype="f")+base.th - - lat=lat.reshape((1,ny,nx)) - lon=lon.reshape((1,ny,nx)) - hgt=hgt.reshape((1,ny,nx)) - + + lat=lat.reshape((1,ny,int(nx))) + lon=lon.reshape((1,ny,int(nx))) + hgt=hgt.reshape((1,ny,int(nx))) + d3dname=("t","z","y","x") ud3dname=("t","z","y","xu") ud2dname=("t","y","xu") @@ -164,13 +164,12 @@ def main(): if fileexists: print("Removing : "+fileexists[0]) os.remove(fileexists[0]) - + io.write(filename, u,varname="U", dims=ud3dname,dtype="f",attributes=dict(units="m/s",description="Horizontal (x) wind speed"), extravars=othervars) if __name__ == '__main__': - global wind_speed, case_study for case in range(3): for ws in [5,10,15,25]: wind_speed=float(ws) diff --git a/helpers/gen_init.py b/helpers/gen_init.py index a078591c..cad72955 100755 --- a/helpers/gen_init.py +++ b/helpers/gen_init.py @@ -26,47 +26,47 @@ def update_base(base,filename,nz): def main(): filename="init" nx,nz,ny=(200.,10.,200) - dims=[nz,ny,nx] - + dims=[int(nz),int(ny),int(nx)] + lonmin=-110.0; lonmax=-100.0; dlon=(lonmax-lonmin)/nx latmin=35.0; latmax=45.0; dlat=(latmax-latmin)/ny base=Bunch(u=10.0,w=0.0,v=0.0, qv=0.0013,qc=0.0, p=100000.0, - th=np.arange(273.0,300,(300-273.0)/nz).reshape((nz,1,1)), + th=np.arange(273.0,300,(300-273.0)/nz).reshape((int(nz),1,1)), dz=400.0) base.z=np.arange(0,nz*base.dz,base.dz) if glob.glob("sounding.txt"): update_base(base,"sounding.txt",nz) nz=base.th.size dims=[nz,ny,nx] - + u=np.zeros(dims,dtype="f")+base.u w=np.zeros(dims,dtype="f")+base.w v=np.zeros(dims,dtype="f")+base.v qv=np.zeros(dims,dtype="f")+base.qv qc=np.zeros(dims,dtype="f")+base.qc coscurve=np.cos(np.arange(dims[2])/dims[2]*2*np.pi+np.pi)+1 - hgt=(coscurve*1000).reshape((1,nx)).repeat(ny,axis=0) - + hgt=(coscurve*1000).reshape((1,int(nx))).repeat(ny,axis=0) + lon=np.arange(lonmin,lonmax,dlon) lat=np.arange(latmin,latmax,dlat) lon,lat=np.meshgrid(lon,lat) - + dz=np.zeros(dims)+base.dz - z=np.zeros(dims,dtype="f")+base.z.reshape((nz,1,1))+hgt.reshape((1,ny,nx)) - + z=np.zeros(dims,dtype="f")+base.z.reshape((int(nz),1,1))+hgt.reshape((1,int(ny),int(nx))) + layer1=(dz[0,:,:]/2) z[0,:,:]+=layer1 for i in range(1,int(nz)): z[i,:,:]=z[i-1,:,:]+(dz[i-1,:,:]+dz[i,:,:])/2.0 - + p=np.zeros(dims,dtype="f")+base.p adjust_p(p,0.0,z) th=np.zeros(dims,dtype="f")+base.th - - + + d3dname=("z","y","x") d2dname=("y","x") othervars=[Bunch(data=v, name="V", dims=d3dname,dtype="f",attributes=dict(units="m/s", description="Horizontal (y) wind speed")), @@ -85,7 +85,7 @@ def main(): if fileexists: print("Removing : "+fileexists[0]) os.remove(fileexists[0]) - + io.write(filename, u,varname="U", dims=d3dname,dtype="f",attributes=dict(units="m/s",description="Horizontal (x) wind speed"), extravars=othervars) diff --git a/helpers/gen_init_ideal.py b/helpers/gen_init_ideal.py index b963eafe..3fa2154e 100755 --- a/helpers/gen_init_ideal.py +++ b/helpers/gen_init_ideal.py @@ -26,7 +26,7 @@ def update_base(base,filename,nz): def main(): filename="init" nx,nz,ny=(20.,10.,19) - dims=[nx,nz,ny] + dims=[int(nx),int(nz),int(ny)] # po=np.log(100000.0) # p1=np.log(50000.0) @@ -55,13 +55,13 @@ def main(): qc=np.zeros(dims,dtype="f")+base.qc coscurve=np.cos(np.arange(dims[0])/dims[0]*2*np.pi+np.pi)+1 # p-=coscurve.reshape((nx,1,1))*15000 - hgt=(coscurve*1000).reshape((nx,1)).repeat(ny,axis=1) + hgt=(coscurve*1000).reshape((int(nx),1)).repeat(ny,axis=1) lon=np.arange(lonmin,lonmax,dlon) lat=np.arange(latmin,latmax,dlat) lat,lon=np.meshgrid(lat,lon) #note that this appears "backwards" but these are C-style, fortran will be reversed dz=np.zeros(dims)+base.dz - z=np.zeros(dims,dtype="f")+base.z.reshape((1,nz,1))+hgt.reshape((nx,1,ny)) + z=np.zeros(dims,dtype="f")+base.z.reshape((1,int(nz),1))+hgt.reshape((int(nx),1,int(ny))) layer1=(dz[:,0,:]/2) z[:,0,:]+=layer1 @@ -70,7 +70,7 @@ def main(): p=np.zeros(dims,dtype="f")+base.p# .reshape((1,nz,1)) adjust_p(p,0.0,z) - th=np.zeros(dims,dtype="f")+base.th.reshape((1,nz,1)) + th=np.zeros(dims,dtype="f")+base.th.reshape((1,int(nz),1)) d3dname=("x","z","y") diff --git a/helpers/wrf/ideal2icar.py b/helpers/wrf/ideal2icar.py index 64be2705..51ad8596 100755 --- a/helpers/wrf/ideal2icar.py +++ b/helpers/wrf/ideal2icar.py @@ -6,7 +6,6 @@ from bunch import Bunch import mygis -import units g=9.81 @@ -34,14 +33,14 @@ def main(inputfile,sounding_file=None): p =mygis.read_nc(inputfile,"P").data.repeat(2,axis=yaxis) + pb phb =mygis.read_nc(inputfile,"PHB").data.repeat(2,axis=yaxis) ph =mygis.read_nc(inputfile,"PH").data.repeat(2,axis=yaxis) + phb - + hgt =mygis.read_nc(inputfile,"HGT").data.repeat(2,axis=yaxis2d) land=mygis.read_nc(inputfile,"XLAND").data.repeat(2,axis=yaxis2d) - + nt,nz,ny,nx=qv.shape print(nx,ny,nz) dims=np.array(qv.shape) - + z=(ph)/g dz=np.diff(z,axis=1) # dz shape = (time,nz-1,ny,nx) @@ -53,7 +52,7 @@ def main(inputfile,sounding_file=None): # wrfz[:,0,...]=dz[:,0,...]/2+hgt # for i in range(1,nz): # wrfz[:,i,:,:]=(dz[:,i,:,:]+dz[:,i-1,:,:])/2+wrfz[:,i-1,:,:] - + mean_dz=dz[0,...].mean(axis=1).mean(axis=1) print("MEAN LEVELS:") print("dz_levels=[") @@ -67,7 +66,7 @@ def main(inputfile,sounding_file=None): else: print("]") - + print("FIRST LEVELS:") print("dz_levels=[") for i in range(0,nz,10): @@ -79,17 +78,17 @@ def main(inputfile,sounding_file=None): print(",".join(curlist)+"]") else: print("]") - - + + # dz=np.zeros(dz.shape)+mean_dz[np.newaxis,:,np.newaxis,np.newaxis] dz=np.zeros(dz.shape)+dz[:,:,:,np.newaxis,0] z=np.zeros(dz.shape) z[:,0,...]=dz[:,0,...]/2+hgt for i in range(1,nz): z[:,i,:,:]=(dz[:,i,:,:]+dz[:,i-1,:,:])/2+z[:,i-1,:,:] - + # adjust_p(p,wrfz,z) - + dx=mygis.read_attr(inputfile,"DX") if type(dx)==np.ndarray: dx=dx[0] @@ -97,7 +96,7 @@ def main(inputfile,sounding_file=None): dlat=dx/111.1 lonmin=-110.0; lonmax=lonmin+nx*dlon latmin=40.0; latmax=latmin+ny*dlat - + udims=copy(dims) udims[-1]+=1 vdims=copy(dims) @@ -114,18 +113,18 @@ def main(inputfile,sounding_file=None): vlon=np.arange(lonmin,lonmax,dlon)[:nx] vlat=np.arange(latmin-dlat/2,latmax+dlat/2,dlat)[:ny+1] vlon,vlat=np.meshgrid(vlon,vlat) - + lat=lat.reshape((1,ny,nx)) lon=lon.reshape((1,ny,nx)) hgt=hgt.reshape((1,ny,nx)) - + d3dname=("t","z","y","x") ud3dname=("t","z","y","xu") ud2dname=("t","y","xu") vd3dname=("t","z","yv","x") vd2dname=("t","yv","x") d2dname=("t","y","x") - + othervars=[Bunch(data=v, name="V", dims=vd3dname,dtype="f",attributes=dict(units="m/s", description="Horizontal (y) wind speed")), # Bunch(data=w, name="W", dims=d3dname, dtype="f",attributes=dict(units="m/s", description="Vertical wind speed")), Bunch(data=qv, name="QVAPOR", dims=d3dname, dtype="f",attributes=dict(units="kg/kg",description="Water vapor mixing ratio")), @@ -152,7 +151,7 @@ def main(inputfile,sounding_file=None): if fileexists: print("Removing : "+fileexists[0]) os.remove(fileexists[0]) - + mygis.write(filename, u,varname="U", dims=ud3dname,dtype="f",attributes=dict(units="m/s",description="Horizontal (x) wind speed"), extravars=othervars) @@ -166,6 +165,5 @@ def main(inputfile,sounding_file=None): sounding_file=sys.argv[2] else: sounding_file=None - - main(filename,sounding_file) + main(filename,sounding_file) diff --git a/helpers/wrf/wind_compare.py b/helpers/wrf/wind_compare.py index 3d392cbd..ae9c9d0d 100755 --- a/helpers/wrf/wind_compare.py +++ b/helpers/wrf/wind_compare.py @@ -5,7 +5,7 @@ import numpy as np import matplotlib.pyplot as plt -from lt_winds import linear_winds +# from lt_winds import linear_winds import ideal_linear import mygis @@ -26,15 +26,15 @@ def load_wrf(filename, preciponly=False): precip/=15.0 # 30 outputsteps = 15 hours if preciponly: return Bunch(precip=precip,hgt=None) - + w=mygis.read_nc(filename,"W").data[time,:,1,:] u=mygis.read_nc(filename,"U").data[time,:,1,:] z=mygis.read_nc(filename,"PH").data[time,:,1,:] z+=mygis.read_nc(filename,"PHB").data[time,:,1,:] z/=9.81 - + hgt=mygis.read_nc(filename,"HGT").data[time,1,:] - + return Bunch(w=w,z=z,hgt=hgt,u=u,precip=precip) def load_icar(filename,h,preciponly=False): @@ -46,35 +46,35 @@ def load_icar(filename,h,preciponly=False): for f in filename: output.append(load_icar(f,h,preciponly=preciponly)) return output - + precip=mygis.read_nc(filename,"rain").data precip=precip[11,1]-precip[10,1] if preciponly: return Bunch(precip=precip) - + u=mygis.read_nc(filename,"u").data[0,1,...] # dudx=-np.diff(u,axis=1) dudx=u[:,:-1]-u[:,1:] - + # w=mygis.read_nc(filename,"w").data[1,...] z=dudx*0 z[0,:]=h+dz_levels[0]/2 nz=z.shape[0] for i in range(1,nz): z[i,:]=z[i-1,:]+(dz_levels[i-1]+dz_levels[i])/2 - + # print(z[:,0]) dz=np.diff(h) # w1=u[:,1:-1]*np.sin(np.arctan(dz[np.newaxis,:]/dx)) # w1=u[:,1:-1]*dz[np.newaxis,:]/np.sqrt(dz[np.newaxis,:]**2 + dx**2) w1=u[:,1:-1]*dz[np.newaxis,:]/dx w1=(w1[:,1:]+w1[:,:-1])/2.0 # * dz_levels[:,np.newaxis]/dx # *10.0 - + w2=dudx*dz_levels[:,np.newaxis]/dx for i in range(1,nz): w2[i,:]+=w2[i-1,:] w2=w2[:,1:-1] - + w=z*0 w[:,1:-1]=w1+w2 # w[:,1:-1]=w2 @@ -82,7 +82,7 @@ def load_icar(filename,h,preciponly=False): def load_linear(zs,T2m=260,u=10.0,v=0,levels=np.array([250,750]),Ndsq=6e-5,dthdz=3.0): """docstring for load_linear""" - + mean_wind=np.round(np.mean(u)) U = mean_wind # dont ask... need to figure this out though # print(U) @@ -91,18 +91,18 @@ def load_linear(zs,T2m=260,u=10.0,v=0,levels=np.array([250,750]),Ndsq=6e-5,dthdz z=0.0 dx,dy=2000.0,2000.0 env_gamma = -dthdz/1000.0 - + if (len(zs)%2) == 1: zs=zs[:-1] - + (Fzs,params) = ideal_linear.get_params(T2m,U,Ndsq,zs,env_gamma) # params.tauf*=2 # print(params) (Pt,w,U3d,z3d)=ideal_linear.solve(Fzs,U,dx,params,zlevels=levels) - + return Bunch(w=w,z=levels,z3d=z3d,hgt=zs,u=U3d,precip=Pt*3600) - - + + # old code for use with lt_winds module, which doesn't seem to work # hgt=hgtin.repeat(2,axis=0) # Ny,Nx=hgt.shape @@ -126,8 +126,8 @@ def load_linear(zs,T2m=260,u=10.0,v=0,levels=np.array([250,750]),Ndsq=6e-5,dthdz # u_out[i,...]=mean_wind+u_hat[0,:] # i+=1 # return Bunch(w=w_out,z=levels,hgt=hgt,u=u_out,precip=p) - - + + def interp2point(w,z,point,verbose=False): """docstring for interp2point""" @@ -150,12 +150,12 @@ def vinterp(d1,d2): for z in range(nz): newdata.w[z,x]=interp2point(w=d1.w[:,x],z=d1.z[:,x],point=d2.z[z,x])#,verbose=(x==nx/2)) # newdata.u[z,x]=interp2point(w=d1.u[:,x],z=d1.z[:,x],point=d2.z[z,x])#,verbose=(x==nx/2)) - + return newdata def main(wrf_file,icar_file,output_file,makeplot=True,plot_legend=False, precip_max=2.0,getPrecip=True,Ndsq=None,t=None,add_file=None,master=False, - text_font_size=8, legend_font_size=8, label_font_size=8, + text_font_size=8, legend_font_size=8, label_font_size=8, draw_xlabels=True, draw_ylabels=True): """docstring for main""" if Ndsq==None: @@ -171,13 +171,13 @@ def main(wrf_file,icar_file,output_file,makeplot=True,plot_legend=False, if add_file: icar2=load_icar(add_file,wrfdata.hgt,preciponly=False) linear_data=load_linear(wrfdata.hgt,T2m=t,u=icardata.u,v=0,levels=icardata.z[:,0],Ndsq=Ndsq,dthdz=3.0) - + if not getPrecip: iwrf=vinterp(wrfdata,icardata) else: iwrf=wrfdata - - + + # lim=10 if makeplot: lim=2 @@ -218,7 +218,7 @@ def main(wrf_file,icar_file,output_file,makeplot=True,plot_legend=False, plt.plot(icardata.w[0,150:250],label="ICAR0",color="blue") # plt.legend(loc=3,ncol=2,fontsize=10) plt.plot([0,100],[0,0],color="black",linestyle="--") - + if getPrecip: if not master: fig=plt.figure(figsize=(6,4.5)) @@ -229,7 +229,7 @@ def main(wrf_file,icar_file,output_file,makeplot=True,plot_legend=False, # precip_max=3.5 # precip_max=2.0 x=np.arange(100)*dx - + plt.plot(x, linear_data.precip[150:250]+offset,color="red",label="Linear",linewidth=precip_width) plt.plot(x, wrfdata.precip[150:250]+offset,color="black",label="WRF",linewidth=precip_width) plt.plot(x, icardata.precip[150:250]+offset,color="green",label="ICAR$_t$",linewidth=precip_width) @@ -251,19 +251,19 @@ def main(wrf_file,icar_file,output_file,makeplot=True,plot_legend=False, label="ICAR$_l$" plt.plot(x, icar2.precip[150:250]+offset,color="blue",label=label,linewidth=precip_width) print("ICAR-l "+str(icar2.precip[150:250].mean())) - + if getPrecip: if master: plot=master else: plot = fig.add_subplot(111) - + case=wrf_file.split("/")[0].split("_")[1:] plt.text(10, (precip_max+offset)*0.75, " U ={0[0]}m/s\n RH={0[1]}\n T ={0[2]}K".format(case),fontsize=text_font_size) # if (case[1]=="0.75") and (case[2]=="260") and (case[0]=="5"): if plot_legend: plt.legend(fontsize=legend_font_size) - + # plot.tick_params(axis='both', which='major', labelsize=0) if draw_xlabels: # if case[0]=="20": @@ -277,26 +277,26 @@ def main(wrf_file,icar_file,output_file,makeplot=True,plot_legend=False, plt.ylabel("Precipitation Rate (mm/hr)", fontsize=label_font_size) else: plot.set_yticklabels("") - + plt.ylim(offset,precip_max+offset) - + # plt.subplot(313) # plt.imshow(wrfdata.w) # plt.colorbar() # plt.title("WRF") # plt.xlim(150,250) - + case=wrf_file.split("/")[0] # print("OUTPUT = wfiles/"+case+"_"+output_file+"_w.png") if not master: plt.savefig("wfiles/"+case+"_"+output_file+"_w.png",dpi=300) plt.close() - + if add_file: return iwrf,icardata, icar2,linear_data else: return iwrf,icardata,linear_data - + bad_cases=[ "wind_15_0.99_280_3", @@ -315,12 +315,12 @@ def mean_precip(Nsquared="1e4"): wrf_file=glob.glob(case+"/wrfout_*")[0] icar_file=glob.glob(case+"/thompson_output_2/icar_out2000_01_01*")[0] wrf,icar,linear=main(wrf_file,icar_file,Nsquared,makeplot=False,getPrecip=True,Ndsq=float(Nsquared.replace("e","e-"))) - + output.append([wrf.precip[150:250].mean(),icar.precip[150:250].mean(),linear.precip[150:250].mean()]) full_data.append([wrf.precip[150:250],icar.precip[150:250]]) else: print("BAD:" + case) - + return np.array(output),full_data if __name__ == '__main__': @@ -330,4 +330,3 @@ def mean_precip(Nsquared="1e4"): main(sys.argv[1],sys.argv[2],sys.argv[3],makeplot=True, getPrecip=False, Ndsq=6e-5,add_file=sys.argv[4]) # main(sys.argv[1],sys.argv[2],sys.argv[3],makeplot=True, getPrecip=False, Ndsq=6e-5,add_file=sys.argv[4]) #main(wrf_file,icar_file,output_file,makeplot=True,getPrecip=True,Ndsq=None,t=None,add_file=None) - \ No newline at end of file diff --git a/helpers/wrf/wrf_vars.py b/helpers/wrf/wrf_vars.py index d3314dec..e248fd73 100644 --- a/helpers/wrf/wrf_vars.py +++ b/helpers/wrf/wrf_vars.py @@ -1,13 +1,13 @@ import operator as op -# note, the elements of wrfvars must either be strings, operators, numbers, or lists. +# note, the elements of wrfvars must either be strings, operators, numbers, or lists. # lists will have their elements loaded one by one # strings will be treated as variable names to load from the file # operators will cause the next element to be loaded operated on with the last element -# numbers will simply be used as is. -# +# numbers will simply be used as is. +# # Example : ["P",op.add,"PB"] means that the variable "P" will be loaded, then the variable "PB" -# the result will be stored in variable "P" in the output file. +# the result will be stored in variable "P" in the output file. steps_per_day=24 # number of output steps per day in the WRF out file @@ -22,7 +22,7 @@ def rename_var(inputvar,dummy): "XLAT_U", "XLONG_U", "XLAT_V", "XLONG_V", ["P",op.add,"PB", rename_var,0], - [["PH",op.add,"PHB"], op.div, 9.8, rename_var, 0], + [["PH",op.add,"PHB"], op.truediv, 9.8, rename_var, 0], "PSFC","HGT", "U","V",["T", op.add, 300, rename_var, 0],"QVAPOR", ["QCLOUD", op.add, "QRAIN"],