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

Last change on this file since 1593 was 1328, checked in by maronga, 11 years ago

last commit documented

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