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

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

last commit documented

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