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

Last change on this file since 1247 was 1110, checked in by maronga, 12 years ago

updates in nc2vdf

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