#!/usr/bin/env python """ Visualization using Mayavi and Matplotlib """ import os, pyproj, Image, sord import numpy as np import matplotlib.pyplot as plt from enthought.mayavi import mlab # parameters title = 'SCEC Chino Hills\nSimulation' format = 'pdf'; dpi = 300.0 format = 'png'; dpi = 150.0 resolution = 'intermediate' resolution = 'high' extent = (-119.8, -115.8), (32.9, 35.1) proj = pyproj.Proj( proj='tmerc', lon_0=-117.8, lat_0=34.0, k=0.001 ) if 1: legend = 'Peak ground velocity' ticklabels = '1', '10', '20 cm/s' ticks = 1, 10, 20 field = 'pv' else: legend = 'Velocity magnitude' ticklabels = '0.5', '10', '20 cm/s' ticks = 0.5, 10, 20 field = 'vm' colormap = [ (0, 0.1, 2, 3), # value (1, 1, 1, 1), # red (0, 0, 1, 1), # green (0, 0, 0, 1), # blue (0, 1, 1, 1), # alpha ] colorexp = 2.0 colorlim = ticks[0], ticks[-1] shape = 125, 250, 800 delta = 1.0, 1.0, 0.1 inches = 6.4, 3.6 ylim = 100.0 xlim = ylim * inches[0] / inches[1] axis = -xlim, xlim, -ylim, ylim pixels = int( dpi * inches[0] ), int( dpi * inches[1] ) point = dpi / 72.0 ppi = 150 # Matplotlib section plt.rc( 'lines', solid_joinstyle='round', solid_capstyle='round', dash_joinstyle='round', dash_capstyle='round', ) plt.close( 0 ) fig0 = plt.figure( 0, inches, dpi=ppi ) ax = fig0.add_axes( [0, 0, 1, 1], axisbg='k', xticks=[], yticks=[] ) for spine in ax.spines.itervalues(): spine.set_color( 'none' ) ax.axis( 'scaled' ) ax.axis( axis ) # cities sites = [ -119.17611, 34.19750, 'baseline', 'center', 'Oxnard', -118.24278, 34.05222, 'baseline', 'center', 'Los Angeles', -118.13583, 34.69806, 'baseline', 'center', 'Lancaster', -117.91361, 33.83528, 'top', 'center', 'Anaheim', -117.37861, 33.19583, 'baseline', 'center', 'Oceanside', -117.28889, 34.10833, 'baseline', 'center', 'San Bernardino', -116.54529, 33.83030, 'baseline', 'center', 'Palm Springs', ] x, y = proj( sites[0::5], sites[1::5] ) z = 20 * np.ones_like( x ) ax.plot( x, y, 'o', ms=2.3, mfc='w', mec='k', mew=1.5, alpha=0.4 ) ax.plot( x, y, 'o', ms=2.3, mfc='w', mec='k', mew=0, alpha=1.0 ) va = sites[2::5] ha = sites[3::5] s = sites[4::5] dy = {'top': -3, 'baseline': 3} for i in range( len( s ) ): sord.plt.text( ax, x[i], y[i]+dy[va[i]], s[i], ha=ha[i], va=va[i], size=7, weight='bold', color='w', edgecolor='k' ) # legend w = 25.0 / (axis[1] - axis[0]) rect = 0.142 - w, 0.08, 2 * w, 0.02 cmap = sord.plt.colormap( colormap, colorexp ) sord.plt.colorbar( fig0, cmap, colorlim, legend, rect, ticks, ticklabels, size=7, weight='bold', color='w', edgecolor='k' ) leg = fig0.add_axes( [0, 0, 1, 1] ) leg.set_axis_off() sord.plt.text( leg, 0.87, 0.95, title, ha='center', va='top', size=10, weight='bold', color='w', edgecolor='k' ) # create overlay if format == 'pdf': over = sord.plt.savefig( fig0, format='pdf', transparent=True, distill=False ) else: aa = 3 mask = sord.plt.savefig( fig0, dpi=aa*dpi, transparent=True ) over = sord.plt.savefig( fig0, dpi=aa*dpi, background='k' ) over[:,:,3] = mask[:,:,3] over = Image.fromarray( over, 'RGBA' ) over = over.resize( pixels, Image.ANTIALIAS ) plt.close( 0 ) # Mayavi section x, y = inches size = int( ppi * x + 2 ), int( ppi * y + 48 ) fig = mlab.figure( 'Viz', size=size ) mlab.clf() fig.scene.disable_render = True fig.scene.set_size( pixels ) fig.scene.render_window.aa_frames = 8 # basemap for kind in 'coastlines', 'borders': x, y = sord.data.mapdata( kind, resolution, extent, 10.0 ) x, y = proj( x, y ) z = 15.0 * np.ones_like( x ) mlab.plot3d( x, y, z, color=(0,0,0), line_width=0.5*point, tube_radius=None ) # topography ddeg = 0.5 / 60.0 z, extent = sord.data.topo( extent ) z *= 0.001 n = z.shape x, y = extent if resolution == 'high': x = x[0] + 0.5 * ddeg * np.arange( n[0] * 2 - 1 ) y = y[0] + 0.5 * ddeg * np.arange( n[1] * 2 - 1 ) z = sord.data.upsample( z ) else: x = x[0] + ddeg * np.arange( n[0] ) y = y[0] + ddeg * np.arange( n[1] ) y, x = np.meshgrid( y, x ) s = np.maximum( 0.01, z ) i = (x + y) < -84.0 s[i] = z[i] x, y = proj( x, y ) cmap = [ (-5, -3, -2, -1, 1, 2, 3, 5), # value ( 0, 0, 10, 10, 15, 15, 25, 25), # red (10, 10, 20, 20, 25, 30, 25, 25), # green (38, 38, 40, 40, 25, 20, 17, 17), # blue (80, 80, 80, 80, 80, 80, 80, 80), # alpha ] cmap = sord.mlab.colormap( cmap, 2.5 ) surf = mlab.mesh( x, y, z, scalars=s, vmin=-4, vmax=4 ) surf.module_manager.scalar_lut_manager.lut.table = cmap surf.actor.property.ambient = 0.0 surf.actor.property.diffuse = 1.0 surf.actor.property.specular = 0.6 surf.actor.property.specular_power = 10 surf.parent.parent.filter.splitting = False # salton sea if 0: ll = (-116.1, -115.5), (33.1, 33.6) x, y = sord.data.mapdata( 'coast', resolution, ll, 800.0, 2, 2 ) x, y = proj( x, y ) x = np.r_[ x.mean(), x ] y = np.r_[ y.mean(), y ] z = np.zeros_like( x ) i = np.arange( 1, x.size-1 ) t = np.array( [np.zeros_like( i ), i, i + 1] ).T surf = mlab.triangular_mesh( x, y, z, t ) surf.module_manager.scalar_lut_manager.lut.table = cmap surf.actor.property.ambient = 0.0 surf.actor.property.diffuse = 1.0 surf.actor.property.specular = 0.6 surf.actor.property.specular_power = 10 surf.parent.parent.filter.splitting = False # lighting fig.scene.light_manager.lights[0].azimuth = 30 fig.scene.light_manager.lights[0].elevation = 30 fig.scene.light_manager.lights[0].intensity = 1.0 fig.scene.light_manager.lights[0].activate = True fig.scene.light_manager.lights[1].activate = False fig.scene.light_manager.lights[2].activate = False # surface plot nn = shape[:2] cmap = sord.mlab.colormap( colormap, colorexp ) x = np.fromfile( 'lon.bin', 'f' ).reshape( nn[::-1] ).T y = np.fromfile( 'lat.bin', 'f' ).reshape( nn[::-1] ).T x, y = proj( x, y ) z = np.zeros_like( x ) surf = mlab.mesh( x, y, z, scalars=z ) surf.module_manager.scalar_lut_manager.lut.table = cmap surf.module_manager.scalar_lut_manager.use_default_range = False surf.module_manager.scalar_lut_manager.data_range = colorlim surf.actor.property.ambient = 0.0 surf.actor.property.diffuse = 0.5 surf.actor.property.specular = 0.3 surf.actor.property.specular_power = 15 surf.parent.parent.filter.splitting = False surf = surf.mlab_source # clock x, y, z, s = ylim * np.array([-1.27, -0.55, 0.2, 0.1]) clock = ( sord.mlab.digital_clock( x, y, z, scale=s, line_width=point*3, color=(0,0,0), opacity=0.2 ), sord.mlab.digital_clock( x, y, z*1.1, scale=s, line_width=point*1.5 ), ) clock[0]( None ) clock[1]( None ) # camera mlab.view( 0, 0, 600, (0, 0, 0) ) fig.scene.parallel_projection = True fig.scene.camera.parallel_scale = axis[3] # combine overlay and save image if field != 'vm': file = 'viz/%s.%s' % (field, format) print file s = np.fromfile( field, 'f' ).reshape( nn[::-1] ).T z = 5 + 5 * s ** 0.1 surf.set( z=z, scalars=s ) fig.scene.disable_render = False out = sord.mlab.screenshot( fig ) if format == 'pdf': out = sord.viz.img2pdf( out, dpi=dpi ) out = sord.viz.pdf_merge( (out, over) ) open( file, 'wb' ).write( out.getvalue() ) else: out = Image.fromarray( out, 'RGB' ) out.paste( over, mask=over ) out.save( file ) else: f1 = open( 'out/v1.bin', 'rb' ) f2 = open( 'out/v2.bin', 'rb' ) f3 = open( 'out/v3.bin', 'rb' ) for i in range( shape[-1] ): file = 'tmp/%s%03d.%s' % (field, i, format) print file x = sord.util.ndread( f1, shape, [(),(),i+1], 'f' ) y = sord.util.ndread( f2, shape, [(),(),i+1], 'f' ) z = sord.util.ndread( f3, shape, [(),(),i+1], 'f' ) s = np.sqrt( x*x + y*y + z*z ) z = 5 + 5 * s ** 0.1 surf.set( z=z, scalars=s ) clock[0]( delta[-1] * i ) clock[1]( delta[-1] * i ) fig.scene.disable_render = False out = sord.mlab.screenshot( fig ) if format == 'pdf': out = sord.viz.img2pdf( out, dpi=dpi ) out = sord.viz.pdf_merge( (out, over) ) open( file, 'wb' ).write( out.getvalue() ) else: out = Image.fromarray( out, 'RGB' ) out.paste( over, mask=over ) out.save( file ) f1.close() f2.close() f3.close()