#! /usr/local/bin/vtkpython import os import sys import pdb from optparse import OptionParser from math import * from vtk_utils import * from netcdf_utils import * from MBARI_utils import * def connect_pipeline(*filts): """ Given a variable number of vtk filters as input, connect them in order. Returns """ lastfilt=filts[0] for iter in filts[1:]: iter.SetInput(lastfilt.GetOutput()) lastfilt=iter return filts def init_pipeline(settings,isovals): """ Sets up the pipeline to use. In order to twiddle around with adding/removing items, we create a bunch of filters, and put the ones we want in a list for easy editing. The 'isovals' argument is a list of isosurface values to generate. This also reads a list of clip planes from the setting argument, instantiates and adds those to the pipeline. """ # Contour filter: the main workhorse of the isosurfacer. con1=vtk.vtkContourFilter() #con1.SetComputeScalars(0) con1.SetComputeNormals(0) # Add the correct isovalues for idx,isoval in enumerate(isovals): con1.SetValue(idx,isoval) # Generate normals from a polygonal mesh. pdn=vtk.vtkPolyDataNormals() pdn.SplittingOff() pdn.ConsistencyOn() #pdn.SetFeatureAngle(90) #pdn.NonManifoldTraversalOn() # Remove duplicate points/remove degenerate cells. clee=vtk.vtkCleanPolyData() clee.PointMergingOff() # Scale depth. scale1=vtk.vtkTransform() scale1.Scale(0.001,1,1) tfm=vtk.vtkTransformPolyDataFilter() tfm.SetTransform(scale1) # Smooth out some of the sharp points. smoo=vtk.vtkSmoothPolyDataFilter() smoo.SetRelaxationFactor(0.1) # Last filter: writes to disk. wtr=vtk.vtkPolyDataWriter() wtr.SetFileTypeToBinary() subd=vtk.vtkButterflySubdivisionFilter() # Converts polygons to triangles. tri=vtk.vtkTriangleFilter() # Converts disjoint triangles to triangle strips for faster rendering stri=vtk.vtkStripper() # Split the pipeline into the filters which come before and after the # clip planes are applied. pipeline_preclip=[con1] pipeline_postclip=[smoo,tfm,tri,pdn,wtr] pipeline=pipeline_preclip[:] # For each clip plane, generate the VTK filter, and append it to the # pipeline. for ((ox,oy,oz),(nx,ny,nz)) in settings.clip_planes: plane=vtk.vtkPlane() plane.SetOrigin(ox,oy,oz) plane.SetNormal(nx,ny,nz) clip=vtk.vtkClipPolyData() clip.SetClipFunction(plane) pipeline.append(clip) pipeline+=pipeline_postclip #String together the filters connect_pipeline(*pipeline) return pipeline def write_iso(sg,vals,pipeline,outfname): # Assume that the last filter in the pipeline is a file writer, and # set the filename. pipeline[-1].SetFileName(outfname) # Attach the pipeline to the data. pipeline[0].SetInput(sg) # Run it. pipeline[-1].Update() if __name__=="__main__": usage="roms_extractor [options] " optparser=OptionParser(usage=usage) optparser.add_option("-u","--u-texture",dest="u_tex", help="Expression to store in the u texture coordinate", type="string", default=None) optparser.add_option("","--u-min",dest="u_min", help="Maximum value for u texture expression", type="float", default=None) optparser.add_option("","--u-max",dest="u_max", help="Minimum value for u texture expression", type="float", default=None) optparser.add_option("-v","--v-texture",dest="v_tex", help="Expression to store in the v texture coordinate", type="string", default=None) optparser.add_option("","--v-min",dest="v_min", help="Maximum value for v texture expression", type="float", default=None) optparser.add_option("","--v-max",dest="v_max", help="Minimum value for v texture expression", type="float", default=None) optparser.add_option("-p","--profile",dest="profile_fname", help="File name of a vertical profile to be used in expressions", type="string", default=None) optparser.add_option("-s","--settings",dest="settings_fname", help="File name of dataset-specific settings", type="string", default="default.settings") optparser.add_option("","--sentinel",dest="sentinel", help="Value to use in place of the sentinel", type="float", default=-9999.0) optparser.add_option("","--mask-var", dest="mask_var", help="Field to use to calculate mask", type="string", default="salinity") optparser.add_option("","--timestep", dest="timestep", help="Index of timestep to use (for 4D datasets)", type="int", default="0") (options,args)=optparser.parse_args() if len(args)!=3: optparser.print_help() sys.exit(2) try: fname=args[0] expr=args[1] isovals=map(float,args[2].split(',')) except: optparser.print_help() sys.exit(2) try: settings=load_settings(options.settings_fname) except Exception,e: print >>sys.stderr,"Error loading settings file:",e sys.exit(0) pipeline=init_pipeline(settings,isovals) ncf=NCFile(fname,settings.name_to_canon) syms={} if (options.profile_fname): syms["mean"],syms["stddev"]=load_profile(options.profile_fname) ncf.set_timestep(options.timestep) # Compute approximate delta_x and delta_y depths=ncf["depth"] lats=ncf["latitude"] longs=ncf["longitude"] # Assume delta(latitude) and delta(longitude) are constant delta_lat=lats[1]-lats[0] delta_long=longs[1]-longs[0] midlat=lats[len(lats)/2] delta_y=earth_circumference(45.0)*delta_lat/360 delta_x=earth_circumference(midlat)*delta_long/360 vel_u=ncf["u"]*settings.velocity_scale vel_v=ncf["v"]*settings.velocity_scale syms["vorticity_z"]=compute_vertical_vorticity(vel_u,vel_v, delta_x,delta_y) missing_value=ncf.missing_value(options.mask_var) mask=Nm.where(ncf[options.mask_var]!=missing_value,1,0) syms["mask"]=mask; isodat=ncf.parse_expression(expr,syms) if options.sentinel: isodat=Nm.where(mask,isodat,options.sentinel) okvals=Nm.compress(Nm.ravel(mask),Nm.ravel(isodat)) print >>sys.stderr,"MIN:",Nm.minimum.reduce(okvals) print >>sys.stderr,"MAX:",Nm.maximum.reduce(okvals) #pdb.set_trace() rg=nm_to_vtk_rectilinear_grid(depths,lats,longs,isodat) u_dat=None v_dat=None if (options.u_tex): u_dat=ncf.parse_expression(options.u_tex,syms) u_dat=normalize_array(u_dat,options.u_min,options.u_max) if (options.v_tex): v_dat=ncf.parse_expression(options.v_tex,syms) v_dat=normalize_array(v_dat,options.v_min,options.v_max) if (u_dat or v_dat): add_textures(rg,u_dat,v_dat) # Run the pipeline with the data, and write it to stdout write_iso(rg,isovals,pipeline,"/dev/fd/1")