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

Last change on this file since 1109 was 1068, checked in by maronga, 12 years ago

last commit documented

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