source: palm/trunk/SCRIPTS/nc2vdf.ncl @ 1046

Last change on this file since 1046 was 1046, checked in by maronga, 9 years ago

put scripts and utilities under GPL

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 13.1 KB
Line 
1;--------------------------------------------------------------------------------!
2; This file is part of PALM.
3;
4; PALM is free software: you can redistribute it and/or modify it under the terms
5; of the GNU General Public License as published by the Free Software Foundation,
6; either version 3 of the License, or (at your option) any later version.
7;
8; PALM is distributed in the hope that it will be useful, but WITHOUT ANY
9; WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
10; A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
11;
12; You should have received a copy of the GNU General Public License along with
13; PALM. If not, see <http://www.gnu.org/licenses/>.
14;
15; Copyright 1997-2012  Leibniz University Hannover
16;------------------------------------------------------------------------------!
17;
18; Current revisions:
19; -----------------
20; Initial revision
21;
22; Former revisions:
23; -----------------
24; $Id: nc2vdf.ncl 1046 2012-11-09 14:38:45Z maronga $
25;
26; Description:
27; ------------
28; This NCL script coverts PALM NetCDF data to vdf, VAPOR's own data format. In
29; order to run this script, NCL version 5.2.0 or higher is required.
30; Default setting will be loaded from .nc2vdf.config.
31; Note that on HLRN, vapor is only available on hicegate2, bicegate2 and the UV
32; sytem.
33;------------------------------------------------------------------------------!
34
35load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
36load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
37load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
38load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
39
40outputfile = "default.vdf"
41
42if(isfilepresent(inputfile)) then
43   print("*** Loading " + inputfile)
44   loadscript(inputfile)
45else
46   print("+++ Input file does not exist. Exiting...")
47   exit
48end if
49
50
51begin
52
53;   print loaded data from input file
54    print("*** loading " + dimsizes(files) + ":")
55    do i=0,dimsizes(files)-1
56       print("    " + i + ") " + files(i)) 
57    end do
58
59    print("*** will create dataset " + outputfile)
60
61
62    workpath = systemfunc("echo $PWD") + "/"
63    path     = ""
64
65
66; hannes stuff below ;-)
67
68;******************************************************
69
70x_min = new(1,double)
71x_max = new(1,double)
72x_min = 99.0
73x_max = -99.0
74y_min = new(1,double)
75y_max = new(1,double)
76y_min = 99.0
77y_max = -99.0
78z_min = new(1,double)
79z_max = new(1,double)
80z_min = 99.0
81z_max = -99.0
82
83vars3d      = ""
84
85
86;print("You may need to run 'vapor-setup.sh' in order to set all enviroment variables.")
87;/muksoft/packages/vapor/2.1.0/bin/vapor-setup.sh
88
89
90;**** Get names of variables wanted to appear in the VDF file, then get their dimensions with min/max x,y,z
91print(" ")
92print("Choose the *.nc files containing 3D data by typing their indices (seperated by ','). Start with the first files of the timeseries.")
93print(" ")
94files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
95print(" "+files_avail)
96
97
98files_in_ary  = toint(str_split(systemfunc("read file; echo $file"), ","))
99nofiles       = dimsizes(files_in_ary)
100file_in       = new(nofiles,string)
101
102do i=0,nofiles-1
103    file_in(i)    = files_avail(files_in_ary(i))
104end do
105
106
107
108
109
110
111if (nofiles .ge. 2) then
112nvars   = 1
113varString  = new((/nofiles,nvars/),string)
114dimNames   = new((/nofiles,nvars,4/),string)
115end if
116
117
118do i=0,nofiles-1
119    f             = addfile(workpath+file_in(i),"r")
120
121    varNames      = getfilevarnames(f)
122    print("Loading file...")
123   
124   
125    print(" ")
126    print("The following variables were found in "+file_in(i)+":")
127    print(" ")
128    print(" "+varNames)
129    print(" ")
130    print("Type the index of the VARIABLES (not dimensions) you want to visualize in Vapor (seperated by ',')")
131
132    vars       = toint(str_split(systemfunc("read vars; echo $vars"), ","))
133
134
135    nvars = dimsizes(vars)
136
137    if (nofiles .eq. 1) then
138    varString  = new((/nofiles,nvars/),string)
139    dimNames   = new((/nofiles,nvars,4/),string)
140    end if
141
142
143
144    do n=0,nvars-1
145
146       varString(i,n) = varNames(vars(n))
147       
148       if (vars3d .ne. "") then
149       vars3d       = vars3d + ":" + varNames(vars(n))
150       else
151       vars3d       = varNames(vars(n))
152       end if
153
154
155       dims          = getfilevardimsizes(f, varString(i,n))
156       dimNames(i,n,:) = getfilevardims(f, varString(i,n))
157
158       
159
160       time     = f->$dimNames(i,n,0)$
161       t_max    = max(time)
162       z        = f->$dimNames(i,n,1)$
163       y        = f->$dimNames(i,n,2)$
164       x        = f->$dimNames(i,n,3)$
165
166       if (z_min .ge. min(z)) then
167          z_min = min(z)
168       end if
169
170       if (z_max .le. max(z)) then
171          z_max = max(z)
172       end if
173
174       if (y_min .ge. min(y)) then
175          y_min = min(y)
176       end if
177
178       if (y_max .le. max(y)) then
179          y_max = max(y)
180       end if
181
182       if (x_min .ge. min(x)) then
183          x_min = min(x)
184       end if
185
186       if (x_max .le. max(x)) then
187          x_max = max(x)
188       end if
189
190       delete(z)
191       delete(y)
192       delete(x)
193    end do
194
195end do
196
197print(varString)
198
199delete(files_avail)
200
201
202print(" ")
203periodic              = systemfunc("read -p 'Periodic lateral boundaries? (y/n) ' RESP; echo $RESP")
204if (periodic .eq. "y") then
205   px           = 1
206   py           = 1
207else
208   px           = 0
209   py           = 0
210end if
211print(" ")
212print(" ")
213
214grid              = systemfunc("read -p 'Is grid stretching used? (y/n) YES will cut of heights above dz_stretch_level? ' RESP; echo $RESP")
215if (grid .eq. "y") then
216   stretch_level    = stringtointeger(systemfunc("read -p 'Please type height of dz_stretch_level (in meters) ' RESP; echo $RESP"))
217   z_max = stretch_level ;**** limit z_max
218end if
219
220print(" ")
221
222;*******************************************************************************
223
224
225
226
227if (any(varNames .EQ. "zwwi")) then
228   var2d     = " -vars2dxy HGT:topo"
229   topo      = True
230   print(" ")
231   print("Found topography - will add 'zwwi' (w-grid) topography to *.vdf file.")
232   print(" ")
233else
234   var2d     = " "
235   topo      = False
236end if
237
238print(" ")
239
240print(dims(0)+" time steps found in "+file_in)
241ts       = systemfunc("read -p 'Include all time steps in *.vdf file? (y/n) ' RESP; echo $RESP")
242
243if (ts .eq. "y") then
244   t_start      = 0
245   t_end        = dimsizes(time)-1
246else
247   t_start      = stringtointeger(systemfunc("read -p 'Type index of first time step: ' ts; echo $ts"))
248   t_end        = stringtointeger(systemfunc("read -p 'Type index of last time step: ' ts; echo $ts"))
249end if
250
251print(" ")
252askreflevel              = systemfunc("read -p 'Specify a maximum refinement level? (y/n) No will set the default value.  ' RESP; echo $RESP")
253if (askreflevel .eq. "y") then
254   reflevel    = systemfunc("read -p 'Please type maximum refinement level ' RESP; echo $RESP")
255   reflevel_str = "-level "+ reflevel
256else
257   reflevel_str = ""
258end if
259print(" ")
260
261
262;*******************************************************************************
263;*******************************************************************************
264;*******************************************************************************
265print(" ")
266print("Creating *.vdf file")
267print(" ")
268;vars3d    = str_join(varString,":")
269
270a = systemfunc(path+"vdfcreate -dimension "+dims(3)+"x"+dims(2)+"x"+dims(1)+" -numts " + (t_end - t_start + 1) + var2d + " "+reflevel_str+" -vars3d "+vars3d+" -extents "+x_min+":"+y_min+":"+z_min+":"+x_max+":"+y_max+":"+z_max+" -periodic "+px+":"+py+":0 "+workpath+outputfile)
271delete(a)
272print("Done.")
273
274
275;*******************************************************************************
276print("Populating *.vdf file")
277
278do i=0,nofiles-1
279    do n=0,nvars-1
280       dimnamesstring = str_join(dimNames(i,n,::-1),":")
281       print(" ")
282       print("Variable: "+varString(i,n)+", number "+(n+1+i)+" of "+(nvars+nofiles-1) )
283       print(" ")
284       do t=t_start,t_end
285          print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
286          a = systemfunc(path+"ncdf2vdf -dimnames "+dimnamesstring+" "+reflevel_str+" -varname "+varString(i,n)+" -ncdfvar "+varString(i,n)+" -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(i))
287          delete(a)
288       end do
289    end do
290end do
291
292if topo then
293   print(" ")
294   print("Adding topography (1)...")
295   print(" ")
296   do t=t_start,t_end
297      print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
298      a = systemfunc(path+"ncdf2vdf -dimnames x:y "+reflevel_str+" -varname HGT -ncdfvar zwwi -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(0))
299      delete(a)
300   end do
301   print(" ")
302   print("Adding topography (2)...")
303   print(" ")
304   do t=t_start,t_end
305      print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
306      a = systemfunc(path+"ncdf2vdf -dimnames x:y "+reflevel_str+" -varname topo -ncdfvar zwwi -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(0))
307      delete(a)
308   end do
309end if
310
311
312;*******************************************************************************
313
314;**** Read more files and add them to VDF/VDC
315addmore = "y"
316do while (addmore .eq. "y")
317
318   addmore       = systemfunc("read -p 'Add time steps or 2D cross sections from additional *.nc files? (2d/time/n) ' RESP; echo $RESP")
319
320;**** Add time steps from additional 3D files
321    if (addmore .eq. "time") then
322
323    delete(f)
324         
325    do i=0,nofiles-1
326
327        print(" ")
328        print("Choose the *.nc file that contains 3D data for "+varString(i,0)+" by typing itÂŽs index.")
329        print(" ")
330        files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
331        print(" "+files_avail)
332        print(" ")
333        file_in(i)    = files_avail(toint(systemfunc("read file; echo $file")))
334       
335    end do
336
337    f        = addfile(workpath+file_in(0),"r")
338
339    time     = f->time
340
341    a = systemfunc(path+"vdfedit -addts "+dimsizes(time)+" "+workpath+outputfile)
342    delete(a)
343
344    do i=0,nofiles-1
345        do n=0,nvars-1
346           dimnamesstring = str_join(dimNames(i,n,::-1),":")
347           print(" ")
348           print("Variable: "+varString(i,n)+", number "+(n+1+i)+" of "+(nvars+nofiles-1) )
349           print(" ")
350           do t=t_end,dimsizes(time)-1
351              print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
352              a = systemfunc(path+"ncdf2vdf -dimnames "+dimnamesstring+" "+reflevel_str+" -varname "+varString(i,n)+" -ncdfvar "+varString(i,n)+" -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(i))
353              delete(a)
354           end do
355        end do
356    end do
357
358    if topo then
359       print(" ")
360       print("Adding topography (1)...")
361       print(" ")
362       do t=t_end,dimsizes(time)-1
363          print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
364          a = systemfunc(path+"ncdf2vdf -dimnames x:y "+reflevel_str+" -varname HGT -ncdfvar zwwi -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(0))
365          delete(a)
366       end do
367       print(" ")
368       print("Adding topography (2)...")
369       print(" ")
370       do t=t_end,dimsizes(time)-1
371          print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
372          a = systemfunc(path+"ncdf2vdf -dimnames x:y "+reflevel_str+" -varname HGT -ncdfvar zwwi -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in(0))
373          delete(a)
374       end do
375    end if
376           
377           
378           
379    t_end       = dimsizes(time)-1
380
381
382   
383    end if
384
385
386;**** Add 2D cross sections
387   if (addmore .eq. "2d") then
388      delete(f)
389      delete(file_in)
390
391      print(" ")
392      print("Choose additional file to add 2D cross sections. Type index.")
393      print(" ")
394      files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
395      print(" "+files_avail)
396
397      file_in       = files_avail(toint(systemfunc("read file; echo $file")))
398      delete(files_avail)
399
400      f             = addfile(workpath+file_in,"r")
401
402      varNames2d      = getfilevarnames(f)
403
404      crosstype       = systemfunc("read -p 'Please specify type of cross section (xy/xz/yz) ' RESP; echo $RESP")
405      print(" ")
406      print("The following variables were found in the file:")
407      print(" ")
408      print(" "+varNames2d)
409      print(" ")
410      print("Type the index of the variables you want to add to the *.vdf file (seperated by ',')")
411
412      vars2d       = toint(str_split(systemfunc("read vars; echo $vars"), ","))
413
414
415
416      nvars2d = dimsizes(vars2d)
417
418      varString2d  = new(dimsizes(vars2d),string)
419      dimNames2d   = new((/dimsizes(vars2d),4/),string)
420
421      do n=0,nvars2d-1
422         varString2d(n) = varNames2d(vars2d(n))
423
424         dims2d     = getfilevardimsizes(f,varString2d(n))
425         dimNames2d(n,:) = getfilevardims(f,varString2d(n))
426         
427
428      end do
429
430      vars2dstr    = str_join(varString2d,":")
431
432      a = systemfunc(path+"vdfedit -addvars2d"+crosstype+" "+vars2dstr+" "+workpath+outputfile)
433      delete(a)
434
435      do n=0,nvars2d-1
436         dimnamesstring = str_join(dimNames2d(n,::-1),":")
437
438         do t=0,dims2d(0)-1
439            a = systemfunc(path+"ncdf2vdf -dimnames "+dimnamesstring+" "+reflevel_str+" -varname "+varNames2d(vars2d(n))+" -ncdfvar "+varNames2d(vars2d(n))+" -cnstnames time -cnstvals "+t+" -ts "+t+" "+workpath+outputfile+" "+workpath+file_in)
440            delete(a)
441         end do
442      end do
443      delete(varNames2d)     
444   end if
445
446   addmore       = systemfunc("read -p 'Add additional *.nc files? (y/n) ' RESP; echo $RESP")
447end do
448
449;*******************************************************************************
450print("Done.")
451;*************
452end
Note: See TracBrowser for help on using the repository browser.