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

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

last commit documented / added nc2vdf

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 13.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;
21; Former revisions:
22; -----------------
23; $Id: nc2vdf.ncl 1047 2012-11-09 15:32:58Z maronga $
24;
25; 1046 2012-11-09 14:38:45Z maronga
26; Initial revision
27;
28; Description:
29; ------------
30; This NCL script coverts PALM NetCDF data to vdf, VAPOR's own data format. In
31; order to run this script, NCL version 5.2.0 or higher is required.
32; Default setting will be loaded from .nc2vdf.config.
33; Note that on HLRN, vapor is only available on hicegate2, bicegate2 and the UV
34; sytem.
35;------------------------------------------------------------------------------!
36
37load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
38load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
39load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
40load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
41
42outputfile = "default.vdf"
43
44if(isfilepresent(inputfile)) then
45   print("*** Loading " + inputfile)
46   loadscript(inputfile)
47else
48   print("+++ Input file does not exist. Exiting...")
49   exit
50end if
51
52
53begin
54
55;   print loaded data from input file
56    print("*** loading " + dimsizes(files) + ":")
57    do i=0,dimsizes(files)-1
58       print("    " + i + ") " + files(i)) 
59    end do
60
61    print("*** will create dataset " + outputfile)
62
63
64    workpath = systemfunc("echo $PWD") + "/"
65    path     = ""
66
67
68; hannes stuff below ;-)
69
70;******************************************************
71
72x_min = new(1,double)
73x_max = new(1,double)
74x_min = 99.0
75x_max = -99.0
76y_min = new(1,double)
77y_max = new(1,double)
78y_min = 99.0
79y_max = -99.0
80z_min = new(1,double)
81z_max = new(1,double)
82z_min = 99.0
83z_max = -99.0
84
85vars3d      = ""
86
87
88;print("You may need to run 'vapor-setup.sh' in order to set all enviroment variables.")
89;/muksoft/packages/vapor/2.1.0/bin/vapor-setup.sh
90
91
92;**** Get names of variables wanted to appear in the VDF file, then get their dimensions with min/max x,y,z
93print(" ")
94print("Choose the *.nc files containing 3D data by typing their indices (seperated by ','). Start with the first files of the timeseries.")
95print(" ")
96files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
97print(" "+files_avail)
98
99
100files_in_ary  = toint(str_split(systemfunc("read file; echo $file"), ","))
101nofiles       = dimsizes(files_in_ary)
102file_in       = new(nofiles,string)
103
104do i=0,nofiles-1
105    file_in(i)    = files_avail(files_in_ary(i))
106end do
107
108
109
110
111
112
113if (nofiles .ge. 2) then
114nvars   = 1
115varString  = new((/nofiles,nvars/),string)
116dimNames   = new((/nofiles,nvars,4/),string)
117end if
118
119
120do i=0,nofiles-1
121    f             = addfile(workpath+file_in(i),"r")
122
123    varNames      = getfilevarnames(f)
124    print("Loading file...")
125   
126   
127    print(" ")
128    print("The following variables were found in "+file_in(i)+":")
129    print(" ")
130    print(" "+varNames)
131    print(" ")
132    print("Type the index of the VARIABLES (not dimensions) you want to visualize in Vapor (seperated by ',')")
133
134    vars       = toint(str_split(systemfunc("read vars; echo $vars"), ","))
135
136
137    nvars = dimsizes(vars)
138
139    if (nofiles .eq. 1) then
140    varString  = new((/nofiles,nvars/),string)
141    dimNames   = new((/nofiles,nvars,4/),string)
142    end if
143
144
145
146    do n=0,nvars-1
147
148       varString(i,n) = varNames(vars(n))
149       
150       if (vars3d .ne. "") then
151       vars3d       = vars3d + ":" + varNames(vars(n))
152       else
153       vars3d       = varNames(vars(n))
154       end if
155
156
157       dims          = getfilevardimsizes(f, varString(i,n))
158       dimNames(i,n,:) = getfilevardims(f, varString(i,n))
159
160       
161
162       time     = f->$dimNames(i,n,0)$
163       t_max    = max(time)
164       z        = f->$dimNames(i,n,1)$
165       y        = f->$dimNames(i,n,2)$
166       x        = f->$dimNames(i,n,3)$
167
168       if (z_min .ge. min(z)) then
169          z_min = min(z)
170       end if
171
172       if (z_max .le. max(z)) then
173          z_max = max(z)
174       end if
175
176       if (y_min .ge. min(y)) then
177          y_min = min(y)
178       end if
179
180       if (y_max .le. max(y)) then
181          y_max = max(y)
182       end if
183
184       if (x_min .ge. min(x)) then
185          x_min = min(x)
186       end if
187
188       if (x_max .le. max(x)) then
189          x_max = max(x)
190       end if
191
192       delete(z)
193       delete(y)
194       delete(x)
195    end do
196
197end do
198
199print(varString)
200
201delete(files_avail)
202
203
204print(" ")
205periodic              = systemfunc("read -p 'Periodic lateral boundaries? (y/n) ' RESP; echo $RESP")
206if (periodic .eq. "y") then
207   px           = 1
208   py           = 1
209else
210   px           = 0
211   py           = 0
212end if
213print(" ")
214print(" ")
215
216grid              = systemfunc("read -p 'Is grid stretching used? (y/n) YES will cut of heights above dz_stretch_level? ' RESP; echo $RESP")
217if (grid .eq. "y") then
218   stretch_level    = stringtointeger(systemfunc("read -p 'Please type height of dz_stretch_level (in meters) ' RESP; echo $RESP"))
219   z_max = stretch_level ;**** limit z_max
220end if
221
222print(" ")
223
224;*******************************************************************************
225
226
227
228
229if (any(varNames .EQ. "zwwi")) then
230   var2d     = " -vars2dxy HGT:topo"
231   topo      = True
232   print(" ")
233   print("Found topography - will add 'zwwi' (w-grid) topography to *.vdf file.")
234   print(" ")
235else
236   var2d     = " "
237   topo      = False
238end if
239
240print(" ")
241
242print(dims(0)+" time steps found in "+file_in)
243ts       = systemfunc("read -p 'Include all time steps in *.vdf file? (y/n) ' RESP; echo $RESP")
244
245if (ts .eq. "y") then
246   t_start      = 0
247   t_end        = dimsizes(time)-1
248else
249   t_start      = stringtointeger(systemfunc("read -p 'Type index of first time step: ' ts; echo $ts"))
250   t_end        = stringtointeger(systemfunc("read -p 'Type index of last time step: ' ts; echo $ts"))
251end if
252
253print(" ")
254askreflevel              = systemfunc("read -p 'Specify a maximum refinement level? (y/n) No will set the default value.  ' RESP; echo $RESP")
255if (askreflevel .eq. "y") then
256   reflevel    = systemfunc("read -p 'Please type maximum refinement level ' RESP; echo $RESP")
257   reflevel_str = "-level "+ reflevel
258else
259   reflevel_str = ""
260end if
261print(" ")
262
263
264;*******************************************************************************
265;*******************************************************************************
266;*******************************************************************************
267print(" ")
268print("Creating *.vdf file")
269print(" ")
270;vars3d    = str_join(varString,":")
271
272a = 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)
273delete(a)
274print("Done.")
275
276
277;*******************************************************************************
278print("Populating *.vdf file")
279
280do i=0,nofiles-1
281    do n=0,nvars-1
282       dimnamesstring = str_join(dimNames(i,n,::-1),":")
283       print(" ")
284       print("Variable: "+varString(i,n)+", number "+(n+1+i)+" of "+(nvars+nofiles-1) )
285       print(" ")
286       do t=t_start,t_end
287          print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
288          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))
289          delete(a)
290       end do
291    end do
292end do
293
294if topo then
295   print(" ")
296   print("Adding topography (1)...")
297   print(" ")
298   do t=t_start,t_end
299      print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
300      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))
301      delete(a)
302   end do
303   print(" ")
304   print("Adding topography (2)...")
305   print(" ")
306   do t=t_start,t_end
307      print("Time step "+(t-t_start+1)+" of "+(t_end-t_start+1) )
308      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))
309      delete(a)
310   end do
311end if
312
313
314;*******************************************************************************
315
316;**** Read more files and add them to VDF/VDC
317addmore = "y"
318do while (addmore .eq. "y")
319
320   addmore       = systemfunc("read -p 'Add time steps or 2D cross sections from additional *.nc files? (2d/time/n) ' RESP; echo $RESP")
321
322;**** Add time steps from additional 3D files
323    if (addmore .eq. "time") then
324
325    delete(f)
326         
327    do i=0,nofiles-1
328
329        print(" ")
330        print("Choose the *.nc file that contains 3D data for "+varString(i,0)+" by typing itÂŽs index.")
331        print(" ")
332        files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
333        print(" "+files_avail)
334        print(" ")
335        file_in(i)    = files_avail(toint(systemfunc("read file; echo $file")))
336       
337    end do
338
339    f        = addfile(workpath+file_in(0),"r")
340
341    time     = f->time
342
343    a = systemfunc(path+"vdfedit -addts "+dimsizes(time)+" "+workpath+outputfile)
344    delete(a)
345
346    do i=0,nofiles-1
347        do n=0,nvars-1
348           dimnamesstring = str_join(dimNames(i,n,::-1),":")
349           print(" ")
350           print("Variable: "+varString(i,n)+", number "+(n+1+i)+" of "+(nvars+nofiles-1) )
351           print(" ")
352           do t=t_end,dimsizes(time)-1
353              print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
354              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))
355              delete(a)
356           end do
357        end do
358    end do
359
360    if topo then
361       print(" ")
362       print("Adding topography (1)...")
363       print(" ")
364       do t=t_end,dimsizes(time)-1
365          print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
366          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))
367          delete(a)
368       end do
369       print(" ")
370       print("Adding topography (2)...")
371       print(" ")
372       do t=t_end,dimsizes(time)-1
373          print("Time step "+(t-t_end+1)+" of "+(dimsizes(time)-t_end) )
374          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))
375          delete(a)
376       end do
377    end if
378           
379           
380           
381    t_end       = dimsizes(time)-1
382
383
384   
385    end if
386
387
388;**** Add 2D cross sections
389   if (addmore .eq. "2d") then
390      delete(f)
391      delete(file_in)
392
393      print(" ")
394      print("Choose additional file to add 2D cross sections. Type index.")
395      print(" ")
396      files_avail   = systemfunc("cd "+workpath+"; ls *.nc")
397      print(" "+files_avail)
398
399      file_in       = files_avail(toint(systemfunc("read file; echo $file")))
400      delete(files_avail)
401
402      f             = addfile(workpath+file_in,"r")
403
404      varNames2d      = getfilevarnames(f)
405
406      crosstype       = systemfunc("read -p 'Please specify type of cross section (xy/xz/yz) ' RESP; echo $RESP")
407      print(" ")
408      print("The following variables were found in the file:")
409      print(" ")
410      print(" "+varNames2d)
411      print(" ")
412      print("Type the index of the variables you want to add to the *.vdf file (seperated by ',')")
413
414      vars2d       = toint(str_split(systemfunc("read vars; echo $vars"), ","))
415
416
417
418      nvars2d = dimsizes(vars2d)
419
420      varString2d  = new(dimsizes(vars2d),string)
421      dimNames2d   = new((/dimsizes(vars2d),4/),string)
422
423      do n=0,nvars2d-1
424         varString2d(n) = varNames2d(vars2d(n))
425
426         dims2d     = getfilevardimsizes(f,varString2d(n))
427         dimNames2d(n,:) = getfilevardims(f,varString2d(n))
428         
429
430      end do
431
432      vars2dstr    = str_join(varString2d,":")
433
434      a = systemfunc(path+"vdfedit -addvars2d"+crosstype+" "+vars2dstr+" "+workpath+outputfile)
435      delete(a)
436
437      do n=0,nvars2d-1
438         dimnamesstring = str_join(dimNames2d(n,::-1),":")
439
440         do t=0,dims2d(0)-1
441            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)
442            delete(a)
443         end do
444      end do
445      delete(varNames2d)     
446   end if
447
448   addmore       = systemfunc("read -p 'Add additional *.nc files? (y/n) ' RESP; echo $RESP")
449end do
450
451;*******************************************************************************
452print("Done.")
453;*************
454end
Note: See TracBrowser for help on using the repository browser.