#!/usr/bin/env python # Script to generate TVTK files from CitcomS hdf files # author: Martin Weier # Copyright (C) 2006 California Institue of Technology # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA #import scipy import sys from datetime import datetime from getopt import getopt, GetoptError from pprint import * from math import * import tables #For HDF support import numpy import pyvtk import sys # defaults path = "./example0.h5" vtk_path = "./vtk_output" vtkfile = "%s.%d.vtk" initial = 0 timesteps= None create_topo = False create_bottom = False create_surface = False create_ascii = False nx = None ny = None nz = None nx_redu=None ny_redu=None nz_redu=None el_nx_redu = None el_ny_redu = None el_nz_redu = None radius_inner = None radius_outer = None nproc_surf = None #Filehandler to the HDF file f = None ##################### polygons3d = [] # arrays containing connectivity information polygons2d = [] counter=0 #Counts iterations of citcom2vtk def print_help(): print "Program to convert CitcomS HDF to Vtk files.\n" print "-p, --path [path to hdf] \n\t Specify input file." print "-o, --output [output filename] \n\t Specify the path to the folder for output files." print ("-i, --initial [initial timestep] \n\t Specify initial timestep to export. If not \n \ \t specified script starts exporting from timestep 0.") print "-t, --timestep [max timestep] \n\t Specify to which timestep you want to export. If not\n \ \t specified export all all timestep starting from intial timestep." print "-x, --nx_reduce [nx] \n\t Set new nx to reduce output grid." print "-y, --ny_reduce [ny] \n\t Set new ny to reduce output grid." print "-z, --nz_reduce [nz] \n\t Set new nz to reduce output grid." print "-b, --bottom \n\t Set to export Bottom information to Vtk file." print "-s, --surface \n\t Set to export Surface information to Vtk file." print "-c, --createtopo \n\t Set to create topography information in bottom and surface Vtk file." print "-a, --ascii \n\t Create Vtk ASCII encoded files instead if binary." print "-h, --help, -? \n\t Print this help." #Iterator for CitcomDataRepresentation(yxz) to VTK(xyz) def vtk_iter(nx,ny,nz): for i in xrange(nx): for j in xrange(ny): for k in xrange(nz): yield k + i * nz + j * nz * nx #Reduces the CitcomS grid def reduce_iter(n,nredu): i=0 n_f=float(n) nredu_f=float(nredu) fl=(n_f-1)/nredu_f redu = 0 for i in xrange(nredu+1): yield int(round(redu)) redu = redu + fl #Transform Vectors in Spherical to Cartesian Coordinates 2d #def velocity2cart2d(vel_colat, vel_lon,x , y): # x1 = vel_colat*cos(x)*cos(y)-vel_lon*sin(y) # y1 = vel_colat*cos(x)*sin(y)+vel_lon*cos(y) # z1 = -vel_colat*sin(x) # return x1,y1,z1 #Converts Spherical to Carthesian Coordinates 2d #def RTF2XYZ2d(vel_colat, vel_lon): # x1 = sin(vel_colat)*cos(vel_lon) # y1 = sin(vel_colat)*sin(vel_lon) # z1 = cos(vel_colat) # return x1,y1,z1 #Transform Vectors in Spherical to Cartesian Coordinates def velocity2cart(vel_colat,vel_long,r, x, y, z): x1 = r*sin(x)*cos(y)+vel_colat*cos(x)*cos(y)-vel_long*sin(y) y1 = r*sin(x)*sin(y)+vel_colat*cos(x)*sin(y)+vel_long*cos(y) z1 = r*cos(x)-vel_colat*sin(x) return x1, y1, z1 #Converts Spherical to Cartesian Coordinates def RTF2XYZ(thet, phi, r): x = r * sin(thet) * cos(phi) y = r * sin(thet) * sin(phi) z = r * cos(thet) return x, y, z #Reads Citcom Files and creates a VTK File def citcom2vtk(t): print "Timestep:",t benchmarkstr = "" #Assign create_bottom and create_surface to bottom and surface #to make them valid in methods namespace bottom = create_bottom surface = create_surface ordered_points = [] #reset Sequences for points ordered_temperature = [] ordered_velocity = [] ordered_visc = [] #Surface and Bottom Points #Initialize empty sequences surf_vec = [] botm_vec = [] surf_topo = [] surf_hflux = [] botm_topo = [] botm_hflux = [] surf_points = [] botm_points = [] for capnr in xrange(nproc_surf): ###Benchmark Point 1 Start## #start = datetime.now() ############################ print "Processing cap",capnr+1,"of",nproc_surf cap = f.root._f_getChild("cap%02d" % capnr) #Information from hdf #This information needs to be read only once hdf_coords = cap.coord[:] hdf_velocity = cap.velocity[t] hdf_temperature = cap.temperature[t] hdf_viscosity = cap.viscosity[t] ###Benchmark Point 1 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6) ###Benchmark Point 2 Start## #start = datetime.now() ############################ #Create Iterator to change data representation nx_redu_iter = reduce_iter(nx,nx_redu) ny_redu_iter = reduce_iter(ny,ny_redu) nz_redu_iter = reduce_iter(nz,nz_redu) #vtk_i = vtk_iter(el_nx_redu,el_ny_redu,el_nz_redu) # read citcom data - zxy (z fastest) for j in xrange(el_ny_redu): j_redu = ny_redu_iter.next() nx_redu_iter = reduce_iter(nx,nx_redu) for i in xrange(el_nx_redu): i_redu = nx_redu_iter.next() nz_redu_iter = reduce_iter(nz,nz_redu) for k in xrange(el_nz_redu): k_redu = nz_redu_iter.next() colat, lon, r = map(float,hdf_coords[i_redu,j_redu,k_redu]) x_coord, y_coord, z_coord = RTF2XYZ(colat,lon,r) ordered_points.append((x_coord,y_coord,z_coord)) ordered_temperature.append(float(hdf_temperature[i_redu,j_redu,k_redu])) ordered_visc.append(float(hdf_viscosity[i_redu,j_redu,k_redu])) vel_colat, vel_lon , vel_r = map(float,hdf_velocity[i_redu,j_redu,k_redu]) x_velo, y_velo, z_velo = velocity2cart(vel_colat,vel_lon,vel_r, colat,lon , r) ordered_velocity.append((x_velo,y_velo,z_velo)) ##Delete Objects for GC del hdf_coords del hdf_velocity del hdf_temperature del hdf_viscosity ###Benchmark Point 2 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6) ###Benchmark Point 3 Start## #start = datetime.now() ############################ #Bottom Information from hdf if bottom == True: try: hdf_bottom_coord = cap.botm.coord[:] hdf_bottom_heatflux = cap.botm.heatflux[t] hdf_bottom_topography = cap.botm.topography[t] hdf_bottom_velocity = cap.botm.velocity[t] except: print "\tCould not find bottom information in file.\n \ Set create bottom to false" bottom = False #Surface Information from hdf if surface==True: try: hdf_surface_coord = cap.surf.coord[:] hdf_surface_heatflux = cap.surf.heatflux[t] hdf_surface_topography = cap.surf.topography[t] hdf_surface_velocity = cap.surf.velocity[t] except: print "\tCould not find surface information in file.\n \ Set create surface to false" surface = False ###Benchmark Point 3 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6) ###Benchmark Point 4 Start## #start = datetime.now() ############################ #Compute surface/bottom topography mean if create_topo: surf_mean=0.0 botm_mean=0.0 if surface: for i in xrange(nx): surf_mean += numpy.mean(hdf_surface_topography[i]) surf_mean = surf_mean/ny if bottom: for i in xrange(nx): botm_mean += numpy.mean(hdf_bottom_topography[i]) botm_mean = botm_mean/nx ###Benchmark Point 4 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6) ###Benchmark Point 5 Start## #start = datetime.now() ############################ #Read Surface and Bottom Data if bottom==True or surface == True: for i in xrange(nx): for j in xrange(ny): if bottom==True: #Bottom Coordinates if create_topo==True: colat, lon = hdf_bottom_coord[i,j] x,y,z = RTF2XYZ(colat,lon,radius_inner+float( (hdf_bottom_topography[i,j]-botm_mean)*(10**21)/(6371000**2/10**(-6))/(3300*10)/1000 )) botm_points.append((x,y,z)) else: colat, lon = hdf_bottom_coord[i,j] x,y,z = RTF2XYZ(colat, lon,radius_inner) botm_points.append((x,y,z)) #Bottom Heatflux botm_hflux.append(float(hdf_bottom_heatflux[i,j])) #Bottom Velocity vel_colat, vel_lon = map(float,hdf_bottom_velocity[i,j]) x,y,z = velocity2cart(vel_colat,vel_lon, radius_inner, colat, lon, radius_inner) botm_vec.append((x,y,z)) if surface==True: #Surface Information if create_topo==True: colat,lon = hdf_surface_coord[i,j] #637100 = Earth radius, 33000 = ? x,y,z = RTF2XYZ(colat,lon,radius_outer+float( (hdf_surface_topography[i,j]-surf_mean)*(10**21)/(6371000**2/10**(-6))/(3300*10)/1000 )) surf_points.append((x,y,z)) else: colat, lon = hdf_surface_coord[i,j] x,y,z = RTF2XYZ(colat, lon,radius_outer) surf_points.append((x,y,z)) #Surface Heatflux surf_hflux.append(float(hdf_surface_heatflux[i,j])) #Surface Velocity vel_colat, vel_lon = map(float,hdf_surface_velocity[i,j]) x,y,z = velocity2cart(vel_colat,vel_lon, radius_outer, colat, lon, radius_outer) surf_vec.append((x,y,z)) #del variables for GC if bottom==True: del hdf_bottom_coord del hdf_bottom_heatflux del hdf_bottom_velocity if surface==True: del hdf_surface_coord del hdf_surface_heatflux del hdf_surface_velocity ###Benchmark Point 5 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6) ###Benchmark Point 6 Start## #start = datetime.now() ############################ ################################################################## #Create Connectivity info if counter==0: #For 3d Data i=1 #Counts X Direction j=1 #Counts Y Direction k=1 #Counts Z Direction for n in xrange((el_nx_redu*el_ny_redu*el_nz_redu)-(el_nz_redu*el_nx_redu)): if (i%el_nz_redu)==0: #X-Values!!! j+=1 #Count Y-Values if (j%el_nx_redu)==0: k+=1 #Count Z-Values if i%el_nz_redu!=0 and j%el_nx_redu!=0: #Check if Box can be created #Get Vertnumbers n0 = n+(capnr*(el_nx_redu*el_ny_redu*el_nz_redu)) n1 = n0+el_nz_redu n2 = n1+el_nz_redu*el_nx_redu n3 = n0+el_nz_redu*el_nx_redu n4 = n0+1 n5 = n4+el_nz_redu n6 = n5+el_nz_redu*el_nx_redu n7 = n4+el_nz_redu*el_nx_redu #Created Polygon Box polygons3d.append([n0,n1,n2,n3,n4,n5,n6,n7]) #Hexahedron VTK Representation i+=1 if bottom==True or surface==True: #Connectivity for 2d-Data i=1 for n in xrange((nx)*(ny) - ny): if i%ny!=0 : n0 = n+(capnr*((nx)*(ny))) n1 = n0+1 n2 = n0+ny n3 = n2+1 polygons2d.append([n0,n1,n2,n3]) i+=1 ###Benchmark Point 6 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf\n" % (delta.seconds + float(delta.microseconds)/1e6) #print benchmarkstr ################################################################# #Write Data to VTK #benchmarkstr = "\n\nIO:\n" ###Benchmark Point 7 Start## #start = datetime.now() ############################ print 'Writing data to vtk...' #Surface Points if surface==True: struct_coords = pyvtk.UnstructuredGrid(surf_points, pixel=polygons2d) #topo_scal = pyvtk.Scalars(surf_topo,'Surface Topography', lookup_table='default') hflux_scal = pyvtk.Scalars(surf_hflux,'Surface Heatflux',lookup_table='default') vel_vec = pyvtk.Vectors(surf_vec,'Surface Velocity Vectors') ## tempdata = pyvtk.PointData(hflux_scal,vel_vec) data = pyvtk.VtkData(struct_coords, tempdata,'CitcomS Output %s Timestep %s' % ('surface info',t)) if create_ascii: data.tofile(vtk_path + (vtkfile % ('surface',t)),) else: data.tofile(vtk_path + (vtkfile % ('surface',t)),'binary') print "Written Surface information to file" ###Benchmark Point 7 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6) ###Benchmark Point 8 Start## #start = datetime.now() ############################ if bottom==True: #Bottom Points struct_coords = pyvtk.UnstructuredGrid(botm_points, pixel=polygons2d) #topo_scal = pyvtk.Scalars(botm_topo,'Bottom Topography','default') hflux_scal = pyvtk.Scalars(botm_hflux,'Bottom Heatflux','default') vel_vec = pyvtk.Vectors(botm_vec,'Bottom Velocity Vectors') ## tempdata = pyvtk.PointData(hflux_scal,vel_vec) data = pyvtk.VtkData(struct_coords, tempdata, 'CitcomS Output %s Timestep %s' % ('Bottom info',t)) if create_ascii: data.tofile(vtk_path + (vtkfile % ('bottom',t))) else: data.tofile(vtk_path + (vtkfile % ('bottom',t)),'binary') print "Written Bottom information to file" ###Benchmark Point 8 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6) ###Benchmark Point 9 Start## #start = datetime.now() #General Data struct_coords = pyvtk.UnstructuredGrid(ordered_points,hexahedron=polygons3d) vel_vec = pyvtk.Vectors(ordered_velocity, 'Velocity Vectors') temp_scal = pyvtk.Scalars(ordered_temperature,'Temperature Scalars','default') visc_scal = pyvtk.Scalars(ordered_visc,'Viscosity Scalars','default') ## tempdata = pyvtk.PointData(temp_scal,visc_scal,vel_vec) data = pyvtk.VtkData(struct_coords, tempdata, 'CitcomS Output %s Timestep:%d NX:%d NY:%d NZ:%d Radius_Inner:%f' % (path,t,el_nx_redu,el_ny_redu,el_nz_redu,radius_inner)) ############################ if create_ascii: data.tofile(vtk_path + (vtkfile % ('general',t))) else: data.tofile(vtk_path + (vtkfile % ('general',t)),'binary') print "Written general data to file" ###Benchmark Point 9 Stop## #delta = datetime.now() - start #benchmarkstr += "%.5lf\n" % (delta.seconds + float(delta.microseconds)/1e6) #print benchmarkstr #print "\n" # parse command line parameters def initialize(): global path global vtk_path global initial global timesteps global create_topo global create_bottom global create_surface global create_ascii global nx global ny global nz global nx_redu global ny_redu global nz_redu global el_nx_redu global el_ny_redu global el_nz_redu global radius_inner global radius_outer global nproc_surf global f try: opts, args = getopt(sys.argv[1:], "p:o:i:t:x:y:z:bscah?", ['path=','output=','timestep=','x=','y=','z=','bottom','surface','createtopo','ascii', 'help','?']) except GetoptError, msg: print "Error: %s" % msg sys.exit(1) if len(opts)<=1: print_help() sys.exit(0) for opt,arg in opts: if opt in ('-p','--path'): path = arg if opt in ('-o','--output'): vtk_path = arg if opt in ('-i','--initial'): try: initial = int(arg) except ValueError: print "Initial is not a number." sys.exit(1) if opt in ('-t','--timestep'): try: timesteps = int(arg) except ValueError: print "Timestep is not a number." sys.exit(1) if opt in ('-x','--nx_reduce'): try: nx_redu = int(arg) except ValueError: print "NX is not a number." if opt in ('-y','--ny_reduce'): try: ny_redu = int(arg) except ValueError: print "NY is not a number." if opt in ('-z','--nz_reduce'): try: nz_redu = int(arg) except ValueError: print "NZ is not a number." if opt in ('-b','--bottom'): create_bottom = True if opt in ('-s','--surface'): create_surface = True if opt in ('-c','--createtopo'): create_topo = True if opt in ('-a','--ascii'): create_ascii = True if opt in ('-h','--help'): print_help() sys.exit(0) if opt == '-?': print_help() sys.exit(0) f = tables.openFile(path,'r') nx = int(f.root.input._v_attrs.nodex) ny = int(f.root.input._v_attrs.nodey) nz = int(f.root.input._v_attrs.nodez) #If not defined as argument read from hdf hdf_timesteps = int(f.root.time.nrows) if timesteps==None or timesteps>hdf_timesteps: timesteps = hdf_timesteps if nx_redu==None: nx_redu = nx-1 if ny_redu==None: ny_redu = ny-1 if nz_redu==None: nz_redu = nz-1 if nx_redu>=nx: nx_redu=nx-1 if ny_redu>=ny: ny_redu=ny-1 if nz_redu>=nz: nz_redu=nz-1 el_nx_redu = nx_redu+1 el_ny_redu = ny_redu+1 el_nz_redu = nz_redu+1 radius_inner = float(f.root.input._v_attrs.radius_inner) radius_outer = float(f.root.input._v_attrs.radius_outer) nproc_surf = int(f.root.input._v_attrs.nproc_surf) ############################################################################### def citcoms_hdf2vtk(): global counter #Call initialize to get and set input params initialize() d1 = datetime.now() print "Converting Hdf to Vtk" print "Initial:",initial, "Timesteps:",timesteps print "NX:",el_nx_redu, "NY:",el_ny_redu, "NZ:", el_nz_redu print "Create Bottom: ",create_bottom, " Create Surface: ", create_surface print "Create Topography: ", create_topo for t in xrange(initial,timesteps): start = datetime.now() citcom2vtk(t) counter+=1 delta = datetime.now() - start print "\t%.3lf sec" % (delta.seconds + float(delta.microseconds)/1e6) d2 = datetime.now() f.close() print "Total: %d seconds" % (d2 - d1).seconds ############################################################################### if __name__ == '__main__': citcoms_hdf2vtk()