-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathflx_ppwrf.py
executable file
·152 lines (131 loc) · 5.83 KB
/
flx_ppwrf.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#!/usr/bin/env python
# This script will accept a set of wrfout files
# as arguments and replace the accumulated precip
# with hourly precipitation data. The input files
# be on the format one file and output time per hour.
#
# Usage: (if put in $PATH)
# ppwrf_hourly_precip.py ./wrfout_d01_*
#
# If ./ contains N files matching 'wrfout_d01_*' then
# the result will be N-1 files named 'pwrfout_d01_*'
# the first of these will have the same output time as
# the second file in the original list.
#-------------------------------------------------------------------------------
from sys import argv # support CLI arguments
#from shutil import copyfile # For creating the output netCDF
from numpy import * # Numeric module
from pycdf import * # NetCDF module
# The pycdf module was chosen due to it's support for
# modifying existing netCDF files which makes the addition
# of a single record variable easier. (We don't need to
# copy the entire file variable by variable)
# pycdf can be installed with support for different numeric modules
# I used numpy because I'm used to it and I have it installed.
# If one wants to use a different numeric module the pycdf package
# has to be reinstalled. (Or one could use both packages but that's
# just messy)
#-------------------------------------------------------------------------------
# The problem:
# WRF outputs accumulated precipiation, we want hourly
# The solution:
# We subtract acc rain of the previous output time
# from the current output time
# Routine for calculating differences between output times and storing this in a
# given variable (ncvar)
def get_precip(vnames,ncvar,f1,f2,fn1,fn2):
for v in vnamesnc:
print(v)
try:
try:
vin1=f1.inq_varid(v)[:] # Load variable data
except CDFError:
print ''.join(['Warning: Variable "',v,'" missing in "',fn_p,
'", skipping'])
raise
try:
vin2=f2.inq_varid(v)[:] # Load variable data
except CDFError:
print ''.join(['Warning: Variable "',v,'" missing in "',fn_c,
'", skipping'])
raise
except CDFError:
pass
continue
#print(vin1)
#print(vin2)
# Add accumulation since last output time to precip_rate
# In this way, the contribution from each precip-type
# is added one-by-one
ncvar[:] = ncvar[:]+(vin2[:]-vin1[:])
#print precip_rate[:]
###MAIN###
# Setup 2 file lists
fn_previous = argv[1:-1] # Previous files
fn_current = argv[2:] # Current files
# Variables we should use for calculating total hourly precipitation:
# If any of the variables are not present in the original files, they will be
# skipped, so there shouldn't be a reason to remove entries.
vnamesnc = ['RAINNC','SNOWNC','GRAUPELNC','HAILNC']
vnamesc = ['RAINC','RAINSH']
for fn_p,fn_c in zip(fn_previous,fn_current): # for each pair of filenames
f1 = CDF(fn_p,NC.NOWRITE)# open last file
f2 = CDF(fn_c,NC.WRITE) # open current file
# Make sure that we are not processing a file which has already been
# processed by this script
try:
test = f2.inq_varid('PRECIP_H') # Attempt to read PRECIP_H
print ''.join(["File ",fn_c," already processed, skipping"])
continue # Move to next iteration (next pair of files)
except CDFError:
pass # If PRECIP_H doesn't exist it is safe to proceed
# Check that we have 1-hour between files
t1 = f1.inq_varid('XTIME')
t2 = f2.inq_varid('XTIME')
#print t2[:]-t1[:]
if t1[:].shape != t2[:].shape:
raise Exception(''.join(["Inconsistent time length of time record in",\
"files: ",fn_p," & ",fn_c,".\n",\
"Both must be of length 1"]) )
elif t2[:].shape[0] != 1:
raise Exception(''.join(["Length of time record in ",fn_c," is not 1"]))
elif t2[:]-t1[:] != 60:
raise Exception(''.join(["Time between ouput files",
fn_p," & ",fn_c,
"is not 60 min"]))
# Set the automode flag so that we don't need to worry about setting the
# define/data mode.
f2.automode()
f2.TITLE = ''.join([f2.TITLE,' POST-PROCESSED TO WORK WITH FLX_WRF2'])
# What we aim to do:
# Check for fields matching our list vnames, any matches
# should be subject to our post processing.
# Predefine hourly precipitation field in current file
precip_rate = f2.def_var('PRECIP_H', NC.FLOAT,
('Time','south_north','west_east') )
precip_rate_conv = f2.def_var('PRECIP_H_CONV', NC.FLOAT,
('Time','south_north','west_east') )
#print precip_rate[:]
# Initialize field as 0
precip_rate[:]=0
precip_rate_conv[:]=0
#print precip_rate[:]
# Set attributes
precip_rate.FieldType = 104
precip_rate.MemoryOrder = 'XY'
precip_rate.description = 'TOTAL LARGE SCALE PRECIPIATION SINCE LAST OUTPUT TIME'
precip_rate.units = 'mm/[OUTPUT TIME INTERVAL]'
precip_rate.stagger = ''
precip_rate.coordinates = 'XLONG XLAT'
precip_rate_conv.FieldType = 104
precip_rate_conv.MemoryOrder = 'XY'
precip_rate_conv.description = 'TOTAL CONVECTIVE PRECIPIATION SINCE LAST OUTPUT TIME'
precip_rate_conv.units = 'mm/[OUTPUT TIME INTERVAL]'
precip_rate_conv.stagger = ''
precip_rate_conv.coordinates = 'XLONG XLAT'
# get hourly precipitation values for convective and nonconvective fields
get_precip(vnamesc,precip_rate,f1,f2,fn_p,fn_c)
get_precip(vnamesnc,precip_rate_conv,f1,f2,fn_p,fn_c)
# Close the files
# (probably not required, but good to clarify that we're done)
del f1,f2