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

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

updates in nc2vdf and associated files

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 16.2 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; Code put under GNU GPL
21;
22; Former revisions:
23; -----------------
24; $Id: nc2vdf.ncl 1067 2012-11-26 10:13:14Z maronga $
25;
26; 1062 2012-11-21 15:57:18Z maronga
27; full batch mode support, loading configuration from an external file
28; fallback to interactive mode if no configuration present
29; batch mode does not support loading more than one variable per file if there
30; is more than one file specified (but there is a workaround)
31; you cannot add time steps to a variable from an additional file in batch mode
32;
33; 1046 2012-11-09 14:38:45Z maronga
34; Initial revision
35;
36; Description:
37; ------------
38; This NCL script coverts PALM NetCDF data to vdf, VAPOR's own data format. In
39; order to run this script, NCL version 5.2.0 or higher is required.
40; Default setting will be loaded from .nc2vdf.config.
41; Note that on HLRN, vapor is only available on hicegate2, bicegate2 and the UV
42; sytem.
43;------------------------------------------------------------------------------!
44
45load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
46load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
47load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
48load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
49
50alt_varnames = False
51man = True
52path = ""
53
54if(isfilepresent(inputfile)) then
55   print("*** Loading " + inputfile)
56   loadscript(inputfile)
57   man = False
58else
59    print("+++ Input file does not exist.")
60    man_resp = systemfunc("read -p 'Enter manual mode? (y/n) ' RESP; echo $RESP")
61    if (man_resp .eq. "y") then
62       man = True
63    else
64       print("+++ Exiting...")
65       exit
66    end if
67end if
68
69
70begin
71
72if (.not. man) then
73;   print loaded data from input file
74    print("*** loading " + dimsizes(files) + ":")
75    do i=0,dimsizes(files)-1
76       print("    " + i + ") " + files(i)) 
77    end do
78
79    print("*** will create dataset " + outputfile)
80end if
81
82;******************************************************
83;give path to the directory where NetCDF files are stored:
84    print("*** Checking hostname...")
85    sys        = "unknown"
86    query      = systemfunc("echo $HOSTNAME")
87    check_muk  = isStrSubset(query,"muk.uni-hannover.de")
88    check_hice = isStrSubset(query,"hicegate2")
89    check_bice = isStrSubset(query,"bicegate2")
90    check_huv  = isStrSubset(query,"huv")
91    check_buv  = isStrSubset(query,"buv")
92    check_WS   = isStrSubset(query,"kookaburra") ;*** insert name of your workstation here
93
94
95
96    if ( check_muk .eq. True ) then
97       print("*** nc2vdf will execute at IMUK")
98       sys  = "IMUK"
99    else
100       if ( check_hice .eq. True .or. check_bice .eq. True ) then
101          print("*** nc2vdf will execute on hicegate2/bicegate2")
102          sys = "HLRN"
103       else
104          if ( check_huv .eq. True .or. check_buv .eq. True ) then
105             print("*** nc2vdf will execute on UV")
106             sys = "HLRN"
107             else
108                if ( check_WS .eq. True) then
109                print("*** nc2vdf will execute on a workstation PC/Mac")
110                sys = "workstation"
111                path     = "/Applications/VAPOR.app/Contents/MacOS/" ;*** costomize path for WS
112          end if
113          end if
114       end if
115    end if
116    if ( sys .eq. "unknown" ) then
117       print("+++ unknown system. Exiting...")
118       exit
119    end if
120
121   
122workpath = systemfunc("echo $PWD") + "/"
123
124
125;******************************************************
126;******************************************************
127if man then
128;specify name of output file:
129outputfile      = "output.vdf" ;<<<<------------- !!!
130end if
131;******************************************************
132;******************************************************
133
134x_min = new(1,double)
135x_max = new(1,double)
136x_min = 99.0
137x_max = -99.0
138y_min = new(1,double)
139y_max = new(1,double)
140y_min = 99.0
141y_max = -99.0
142z_min = new(1,double)
143z_max = new(1,double)
144z_min = 99.0
145z_max = -99.0
146
147vars3d      = ""
148
149
150
151if man then
152    ;**** Get names of variables wanted to appear in the VDF file, then get their dimensions with min/max x,y,z
153    print(" ")
154    print("Choose the *.nc files containing 3D data by typing their indices (seperated by ','). Start with the first files of the timeseries.")
155    print(" ")
156    files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
157    print(" "+files_avail)
158
159
160    files_in_ary  = toint(str_split(systemfunc("read file; echo $file"), ","))
161    nofiles       = dimsizes(files_in_ary)
162    file_in       = new(nofiles,string)
163
164    do i=0,nofiles-1
165        file_in(i)    = files_avail(files_in_ary(i))
166    end do
167else
168    file_in = files
169    nofiles = dimsizes(file_in)
170end if
171
172
173
174
175
176if (nofiles .ge. 2) then
177nvars   = 1
178if man then
179varString  = new((/nofiles,nvars/),string)
180end if; man
181dimNames   = new((/nofiles,nvars,4/),string)
182end if
183
184
185do i=0,nofiles-1
186    print("Loading file...")
187    f             = addfile(workpath+file_in(i),"r")
188
189    varNames      = getfilevarnames(f)
190   
191   
192    if man then
193    print(" ")
194    print("The following variables were found in "+file_in(i)+":")
195    print(" ")
196    print(" "+varNames)
197    print(" ")
198    print("Type the index of the VARIABLES (not dimensions) you want to visualize in Vapor (seperated by ',')")
199
200    vars       = toint(str_split(systemfunc("read vars; echo $vars"), ","))
201
202
203    nvars = dimsizes(vars)
204    end if; man
205   
206    if (nofiles .eq. 1) then
207    if man then
208    varString  = new((/nofiles,nvars/),string)
209    end if; man
210    dimNames   = new((/nofiles,nvars,4/),string)
211    end if
212
213
214
215    do n=0,nvars-1
216       if man then
217       varString(i,n) = varNames(vars(n))
218       end if; man
219       
220       if alt_varnames then
221            varString(i,n) = varNames(vars(i,n))
222       end if
223       
224       if (vars3d .ne. "") then
225       vars3d       = vars3d + ":" + varString(i,n)
226       else
227       vars3d       = varString(i,n)
228       end if
229
230
231       dims          = getfilevardimsizes(f, varString(i,n))
232       dimNames(i,n,:) = getfilevardims(f, varString(i,n))
233
234       
235
236       time     = f->$dimNames(i,n,0)$
237       t_max    = max(time)
238       z        = f->$dimNames(i,n,1)$
239       y        = f->$dimNames(i,n,2)$
240       x            = f->$dimNames(i,n,3)$
241
242       if (z_min .ge. min(z)) then
243          z_min = min(z)
244       end if
245
246       if (z_max .le. max(z)) then
247          z_max = max(z)
248       end if
249
250       if (y_min .ge. min(y)) then
251          y_min = min(y)
252       end if
253
254       if (y_max .le. max(y)) then
255          y_max = max(y)
256       end if
257
258       if (x_min .ge. min(x)) then
259          x_min = min(x)
260       end if
261
262       if (x_max .le. max(x)) then
263          x_max = max(x)
264       end if
265
266       delete(z)
267       delete(y)
268       delete(x)
269    end do
270
271end do
272
273;print(varString)
274
275if man then
276delete(files_avail)
277end if
278
279if man then
280    print(" ")
281    periodic_resp              = systemfunc("read -p 'Periodic lateral boundaries? (y/n) ' RESP; echo $RESP")
282    if (periodic_resp .eq. "y") then
283        periodic = True
284    else
285        periodic = False
286    end if
287end if
288
289if (periodic) then
290   px           = 1
291   py           = 1
292else
293   px           = 0
294   py           = 0
295end if
296print(" ")
297print(" ")
298
299if man then
300    grid_resp              = systemfunc("read -p 'Is grid stretching used? (y/n) YES will cut of heights above dz_stretch_level? ' RESP; echo $RESP")
301    if (grid_resp .eq. "y") then
302        grid_stretch = True
303    else
304        grid_stretch = False
305    end if
306end if
307
308if (grid_stretch) then
309    if man then
310       stretch_level    = stringtointeger(systemfunc("read -p 'Please type height of dz_stretch_level (in meters) ' RESP; echo $RESP"))
311    end if; man
312    z_max = stretch_level ;**** limit z_max
313end if
314
315print(" ")
316
317;*******************************************************************************
318
319
320
321
322if (any(varNames .EQ. "zwwi")) then
323   var2d     = " -vars2dxy HGT:topo"
324   topo      = True
325   print(" ")
326   print("Found topography - will add 'zwwi' (w-grid) topography to *.vdf file.")
327   print(" ")
328else
329   var2d     = " "
330   topo      = False
331end if
332
333if man then
334print(" ")
335
336print(dims(0)+" time steps found in "+file_in)
337ts       = systemfunc("read -p 'Include all time steps in *.vdf file? (y/n) ' RESP; echo $RESP")
338
339if (ts .eq. "y") then
340   t_start      = 0
341   t_end        = dimsizes(time)-1
342else
343   t_start      = stringtointeger(systemfunc("read -p 'Type index of first time step: ' ts; echo $ts"))
344   t_end        = stringtointeger(systemfunc("read -p 'Type index of last time step: ' ts; echo $ts"))
345end if
346end if
347
348if (.not. man) then
349    if (.not. spec_ts) then
350           t_start      = 0
351           t_end        = dimsizes(time)-1
352    end if
353end if
354
355if man then
356    print(" ")
357    askreflevel              = systemfunc("read -p 'Specify a maximum refinement level? (y/n) No will set the default value.  ' RESP; echo $RESP")
358    if (askreflevel .eq. "y") then
359       reflevel    = systemfunc("read -p 'Please type maximum refinement level ' RESP; echo $RESP")
360       reflevel_str = "-level "+ reflevel
361    else
362       reflevel_str = ""
363    end if
364    print(" ")
365else
366    if spec_reflevel then
367        reflevel_str = "-level "+ reflevel
368    else
369       reflevel_str = ""
370    end if
371end if
372
373       
374
375;*******************************************************************************
376;*******************************************************************************
377;*******************************************************************************
378print(" ")
379print("Creating *.vdf file")
380print(" ")
381;vars3d    = str_join(varString,":")
382
383a = 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)
384delete(a)
385print("Done.")
386
387
388;*******************************************************************************
389print("Populating *.vdf file")
390
391do i=0,nofiles-1
392    do n=0,nvars-1
393       dimnamesstring = str_join(dimNames(i,n,::-1),":")
394       print(" ")
395       print("Variable: "+varString(i,n)+", number "+(n+1+i)+" of "+(nvars+nofiles-1) )
396       print(" ")
397       do t=t_start,t_end
398          print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
399          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))
400          delete(a)
401       end do
402    end do
403end do
404
405if topo then
406   print(" ")
407   print("Adding topography (1)...")
408   print(" ")
409   do t=t_start,t_end
410      print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
411      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))
412      delete(a)
413   end do
414   print(" ")
415   print("Adding topography (2)...")
416   print(" ")
417   do t=t_start,t_end
418      print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
419      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))
420      delete(a)
421   end do
422end if
423
424
425;*******************************************************************************
426if man then
427;**** Read more files and add them to VDF/VDC
428addmore = "y"
429do while (addmore .eq. "y")
430
431   addmore       = systemfunc("read -p 'Add time steps or 2D cross sections from additional *.nc files? (2d/time/n) ' RESP; echo $RESP")
432
433;**** Add time steps from additional 3D files
434    if (addmore .eq. "time") then
435
436    delete(f)
437         
438    do i=0,nofiles-1
439
440        print(" ")
441        print("Choose the *.nc file that contains 3D data for "+varString(i,0)+" by typing itÂŽs index.")
442        print(" ")
443        files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
444        print(" "+files_avail)
445        print(" ")
446        file_in(i)    = files_avail(toint(systemfunc("read file; echo $file")))
447       
448    end do
449
450    f        = addfile(workpath+file_in(0),"r")
451
452    time     = f->time
453
454    a = systemfunc(path+"vdfedit -addts "+dimsizes(time)+" "+workpath+outputfile)
455    delete(a)
456
457    do i=0,nofiles-1
458        do n=0,nvars-1
459           dimnamesstring = str_join(dimNames(i,n,::-1),":")
460           print(" ")
461           print("Variable: "+varString(i,n)+", number "+(n+1+i)+" of "+(nvars+nofiles-1) )
462           print(" ")
463           do t=t_end,dimsizes(time)-1
464              print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
465              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))
466              delete(a)
467           end do
468        end do
469    end do
470
471    if topo then
472       print(" ")
473       print("Adding topography (1)...")
474       print(" ")
475       do t=t_end,dimsizes(time)-1
476          print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
477          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))
478          delete(a)
479       end do
480       print(" ")
481       print("Adding topography (2)...")
482       print(" ")
483       do t=t_end,dimsizes(time)-1
484          print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
485          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))
486          delete(a)
487       end do
488    end if
489           
490           
491           
492    t_end       = dimsizes(time)-1
493
494
495   
496    end if
497
498
499;**** Add 2D cross sections
500   if (addmore .eq. "2d") then
501      delete(f)
502      delete(file_in)
503
504      print(" ")
505      print("Choose additional file to add 2D cross sections. Type index.")
506      print(" ")
507      files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
508      print(" "+files_avail)
509
510      file_in       = files_avail(toint(systemfunc("read file; echo $file")))
511      delete(files_avail)
512
513      f             = addfile(workpath+file_in,"r")
514
515      varNames2d      = getfilevarnames(f)
516
517      crosstype       = systemfunc("read -p 'Please specify type of cross section (xy/xz/yz) ' RESP; echo $RESP")
518      print(" ")
519      print("The following variables were found in the file:")
520      print(" ")
521      print(" "+varNames2d)
522      print(" ")
523      print("Type the index of the variables you want to add to the *.vdf file (seperated by ',')")
524
525      vars2d       = toint(str_split(systemfunc("read vars; echo $vars"), ","))
526
527
528
529      nvars2d = dimsizes(vars2d)
530
531      varString2d  = new(dimsizes(vars2d),string)
532      dimNames2d   = new((/dimsizes(vars2d),4/),string)
533
534      do n=0,nvars2d-1
535         varString2d(n) = varNames2d(vars2d(n))
536
537         dims2d     = getfilevardimsizes(f,varString2d(n))
538         dimNames2d(n,:) = getfilevardims(f,varString2d(n))
539         
540
541      end do
542
543      vars2dstr    = str_join(varString2d,":")
544
545      a = systemfunc(path+"vdfedit -addvars2d"+crosstype+" "+vars2dstr+" "+workpath+outputfile)
546      delete(a)
547
548      do n=0,nvars2d-1
549         dimnamesstring = str_join(dimNames2d(n,::-1),":")
550
551         do t=0,dims2d(0)-1
552            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)
553            delete(a)
554         end do
555      end do
556      delete(varNames2d)     
557   end if
558
559   addmore       = systemfunc("read -p 'Add additional *.nc files? (y/n) ' RESP; echo $RESP")
560end do
561end if; man
562;*******************************************************************************
563print("Done.")
564;*************
565end
Note: See TracBrowser for help on using the repository browser.