|
| 1 | +from pylab import * |
| 2 | +import glob, os |
| 3 | +import pyvista as pv |
| 4 | +from clawpack.geoclaw import topotools |
| 5 | +from clawpack.visclaw import animation_tools |
| 6 | +from clawpack.geoclaw import fgout_tools |
| 7 | +from clawpack.visclaw import colormaps |
| 8 | + |
| 9 | +warpfactor = 5 # vertical amplification of elevations in 3D plots |
| 10 | + |
| 11 | +outdir = '_output' |
| 12 | +fgno = 1 # fgout grid number |
| 13 | + |
| 14 | +if 1: |
| 15 | + # determine how many fgout frames are in outdir: |
| 16 | + fgout_frames = glob.glob(os.path.join(outdir, \ |
| 17 | + 'fgout%s.t*' % str(fgno).zfill(4))) |
| 18 | + nfgout = len(fgout_frames) |
| 19 | + fgframenos = array(range(1, nfgout+1)) |
| 20 | + print('Found %i fgout frames' % nfgout) |
| 21 | +else: |
| 22 | + # specify which frames to use: |
| 23 | + fgframenos = array(range(1,50,2)) |
| 24 | + |
| 25 | + |
| 26 | +# where to save a png file for each frame, for constructing an animation: |
| 27 | +framedir = '_frames' |
| 28 | +os.system('mkdir -p %s' % framedir) |
| 29 | +os.system('rm %s/*' % framedir) # remove frames from previous version |
| 30 | + |
| 31 | +window_size = (1200,1200) |
| 32 | + |
| 33 | +tan = [0.8,0.5,0.2] |
| 34 | +light_blue = [0.0,1.0,1.0] |
| 35 | +cmap_massfrac = colormaps.make_colormap({0:light_blue, 0.1:light_blue, |
| 36 | + 1:tan}) |
| 37 | + |
| 38 | +def plot_topo(): |
| 39 | + """ |
| 40 | + plot topo from topofile |
| 41 | + """ |
| 42 | + topo = topotools.Topography('basal_topo.tt3') |
| 43 | + z = array([0.]) |
| 44 | + x = topo.x |
| 45 | + y = -topo.y |
| 46 | + X,Y,Z = meshgrid(x, y, z, indexing='ij') |
| 47 | + topoxyz = pv.StructuredGrid(X,Y,Z) |
| 48 | + |
| 49 | + B = flipud(topo.Z) |
| 50 | + B = fliplr(B) |
| 51 | + topoxyz.point_data['B'] = B.flatten(order='C') |
| 52 | + |
| 53 | + topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor) |
| 54 | + |
| 55 | + p = pv.Plotter() |
| 56 | + p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400)) |
| 57 | + p.add_title('Topography') |
| 58 | + p.show(window_size=window_size) |
| 59 | + |
| 60 | +def plot_topo_fgout(): |
| 61 | + """ |
| 62 | + plot topo from an fgout file |
| 63 | + """ |
| 64 | + fgno = 1 # which fgout grid |
| 65 | + |
| 66 | + fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format, qmap='dclaw') |
| 67 | + |
| 68 | + fgout = fgout_grid.read_frame(1) |
| 69 | + |
| 70 | + x = fgout.x |
| 71 | + y = fgout.y |
| 72 | + z = array([0.]) |
| 73 | + X,Y,Z = meshgrid(x, y, z, indexing='ij') |
| 74 | + topoxyz = pv.StructuredGrid(X,Y,Z) |
| 75 | + |
| 76 | + B = fgout.B |
| 77 | + topoxyz.point_data['B'] = B.flatten(order='C') |
| 78 | + warpfactor = 5 # amplification of elevations |
| 79 | + topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor) |
| 80 | + |
| 81 | + h = fgout.h |
| 82 | + topoxyz.point_data['h'] = h.flatten(order='C') |
| 83 | + warpfactor = 5 # amplification of elevations |
| 84 | + hwarp = topoxyz.warp_by_scalar('h', factor=warpfactor) |
| 85 | + |
| 86 | + eta = fgout.eta |
| 87 | + topoxyz.point_data['eta'] = B.flatten(order='C') |
| 88 | + warpfactor = 5 # amplification of elevations |
| 89 | + etawarp = topoxyz.warp_by_scalar('eta', factor=warpfactor) |
| 90 | + |
| 91 | + p = pv.Plotter() |
| 92 | + p.window_size = window_size |
| 93 | + p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400)) |
| 94 | + #p.add_mesh(topowarp,color='g') |
| 95 | + p.add_title('Topography fgout.B') |
| 96 | + p.show() |
| 97 | + |
| 98 | + |
| 99 | +def plot_fgout(make_animation=False): |
| 100 | + """ |
| 101 | + If make_animation == False, plot one fgout frame with a slider bar to |
| 102 | + change frame number. Also prints the camera position after changing frame, |
| 103 | + which you can copy and paste into this function in order to set the |
| 104 | + camera position for making an animation. |
| 105 | + |
| 106 | + If make_animation == True, make an mp4 animation of all frames |
| 107 | + specified by fgframenos set above. |
| 108 | + """ |
| 109 | + |
| 110 | + global etamesh |
| 111 | + |
| 112 | + fgno = 1 # which fgout grid |
| 113 | + warpfactor = 5 # amplification of elevations |
| 114 | + |
| 115 | + # Instantiate object for reading fgout frames: |
| 116 | + fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format, qmap='dclaw') |
| 117 | + |
| 118 | + fgout = fgout_grid.read_frame(1) |
| 119 | + |
| 120 | + x = fgout.x |
| 121 | + y = fgout.y |
| 122 | + z = array([0.]) |
| 123 | + X,Y,Z = meshgrid(x, y, z, indexing='ij') |
| 124 | + topoxyz = pv.StructuredGrid(X,Y,Z) |
| 125 | + |
| 126 | + B = fgout.B |
| 127 | + topoxyz.point_data['B'] = B.flatten(order='C') |
| 128 | + topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor) |
| 129 | + |
| 130 | + p = pv.Plotter(off_screen=make_animation, lighting='three lights') |
| 131 | + p.window_size = window_size |
| 132 | + #p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400)) |
| 133 | + p.add_mesh(topowarp,color='lightgreen') |
| 134 | + etamesh = None |
| 135 | + |
| 136 | + |
| 137 | + def set_frameno(fgframeno): |
| 138 | + global etamesh |
| 139 | + fgframeno = int(round(fgframeno)) |
| 140 | + fgout = fgout_grid.read_frame(fgframeno) |
| 141 | + tsec = fgout.t |
| 142 | + print('Frame %i, t = %.1f seconds' % (fgframeno, fgout.t)) |
| 143 | + |
| 144 | + # replace land surface by nan in eta so it only shows water: |
| 145 | + # (big tolerance for landslide, would normally be smaller) |
| 146 | + eta = where(fgout.h>1, fgout.eta, nan) |
| 147 | + |
| 148 | + topoxyz.point_data['eta'] = eta.flatten(order='C') |
| 149 | + etawarp = topoxyz.warp_by_scalar('eta', factor=warpfactor) |
| 150 | + if etamesh: |
| 151 | + p.remove_actor(etamesh) |
| 152 | + |
| 153 | + # color the mesh based on mass fraction, using cmap_massfrac |
| 154 | + h = fgout.h |
| 155 | + massfrac = divide(fgout.hm, h, where=h>0, out=nan*ones(h.shape)) |
| 156 | + etamesh = p.add_mesh(etawarp,scalars=massfrac, |
| 157 | + colormap=cmap_massfrac, clim=(0,0.6)) |
| 158 | + |
| 159 | + p.add_title('Frame %i at time %.1f seconds' % (fgframeno,tsec)) |
| 160 | + |
| 161 | + if not make_animation: |
| 162 | + # print camera position so that this can be copied and pasted |
| 163 | + # into this script after adjusting (and then sliding frameno) |
| 164 | + print('p.camera_position = ', p.camera_position) |
| 165 | + |
| 166 | + |
| 167 | + # initial camera position: |
| 168 | + p.camera_position = [ |
| 169 | + (9034.796206317897, -2484.6171987620633, 3703.6128928929047), |
| 170 | + (1500.0, 1500.0, 950.89111328125), |
| 171 | + (-0.29984440650740857, 0.08916816605502331, 0.949811755059182)] |
| 172 | + |
| 173 | + |
| 174 | + if not make_animation: |
| 175 | + fgfr1 = fgframenos[0] |
| 176 | + fgfr2 = fgframenos[-1] |
| 177 | + p.add_slider_widget(set_frameno, [fgfr1,fgfr2], |
| 178 | + value=fgfr1, title='Frame', |
| 179 | + pointa=(0.4,0.85), pointb=(0.9,0.85), color='blue', |
| 180 | + slider_width=0.02, tube_width=0.005) |
| 181 | + |
| 182 | + p.show() |
| 183 | + |
| 184 | + else: |
| 185 | + |
| 186 | + # make a png file for each frame: |
| 187 | + for fgframeno in fgframenos: |
| 188 | + set_frameno(fgframeno) |
| 189 | + fname_png = '%s/PyVistaFrame%s.png' \ |
| 190 | + % (framedir, str(fgframeno).zfill(4)) |
| 191 | + p.screenshot(fname_png) |
| 192 | + print('Created ',fname_png) |
| 193 | + |
| 194 | + p.close() |
| 195 | + |
| 196 | + # combine png files into mp4 animation: |
| 197 | + anim = animation_tools.make_anim(framedir, |
| 198 | + fname_pattern='PyVistaFrame*.png') |
| 199 | + fname_mp4 = 'radial_slide_animation.mp4' |
| 200 | + animation_tools.make_mp4(anim, fname_mp4) |
| 201 | + print('Created ',fname_mp4) |
| 202 | + |
| 203 | +if __name__=='__main__': |
| 204 | + |
| 205 | + plot_fgout(make_animation=False) |
0 commit comments