From c8f627a8e3dad20562a49f1662432a3d0b617b07 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. --- helpers/README.md | 19 ++++++++-- 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, 148 insertions(+), 142 deletions(-) diff --git a/helpers/README.md b/helpers/README.md index 238c73c1..d5137e73 100644 --- a/helpers/README.md +++ b/helpers/README.md @@ -1,15 +1,28 @@ # 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 +$ conda env config vars set PYTHONPATH=$(pwd)/lib:$PYTHONPATH + +reactivate environment, this will be saved for future use +$ 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. 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"],