plot_ocean.py.html | ![]() |
Source file: plot_ocean.py | |
Directory: /Users/rjl/clawpack_src/clawpack_master/geoclaw/examples/tsunami/radial-ocean-island-fgmax | |
Converted: Mon Feb 10 2020 at 10:51:30 using clawcode2html | |
This documentation file will not reflect any later changes in the source file. |
from __future__ import print_function from pylab import * from clawpack.geoclaw import topotools from clawpack.geoclaw.data import Rearth # radius of earth from clawpack.clawutil.data import ClawData from mapper import latlong import maketopo probdata = ClawData() probdata.read('setprob.data', force=True) theta_island = probdata.theta_island print("theta_island = ",theta_island) #close(1) #figure(1, figsize=(10,6)) #axes([.1,.1,.5,.8]) clf() s = linspace(0, 360., 1001) xshore,yshore = latlong(1645.e3, s, 40., Rearth) plot(xshore,yshore,'k',linewidth=2) xshelf,yshelf = latlong(1560.e3, s, 40., Rearth) plot(xshelf,yshelf,'k--',linewidth=2) theta_island = 220. (xi,yi) = latlong(1600.e3, theta_island, 40., Rearth) xb = [xi-1, xi+1, xi+1, xi-1, xi-1] yb = [yi-1, yi-1, yi+1, yi+1, yi-1] plot(xb,yb,'b',linewidth=1.5) text(5,49,'Test 1', fontsize=20) theta_island = 260. (xi,yi) = latlong(1600.e3, theta_island, 40., Rearth) xb = [xi-1, xi+1, xi+1, xi-1, xi-1] yb = [yi-1, yi-1, yi+1, yi+1, yi-1] plot(xb,yb,'b',linewidth=1.5) text(9,40,'Test 2', fontsize=20) title('(a) Radial ocean',fontsize=20) x = linspace(-5,5,51) y = linspace(35,45,51) X,Y = meshgrid(x,y) Z = maketopo.qinit(X,Y) contour(X,Y,Z,[2],linewidth=2,colors='k') axis('scaled') xlim([-21,21]) ylim([20,60]) xticks(fontsize='15') yticks([25,35,45,55],fontsize='15') xlabel('Longitude',fontsize=15) ylabel('Latitude',fontsize=15) savefig('ocean.png') print("Created ocean.png") if 0: axes([.65,.6,.25,.25]) plot(xshore,yshore,'k',linewidth=2) plot(xshelf,yshelf,'k--',linewidth=2) maketopo.theta_island = 220. (xisland,yisland) = latlong(1600.e3, maketopo.theta_island, 40., Rearth) x = linspace(xisland-1, xisland+1, 200) y = linspace(yisland-1, yisland+1, 200) X,Y = meshgrid(x,y) Z = maketopo.topo(X,Y) contour(X,Y,Z,[0,7,15]) axis('scaled') xlim([xisland-0.5, xisland+0.5]) ylim([yisland-0.5, yisland+0.5]) axes([.65,.2,.25,.25]) plot(xshore,yshore,'k',linewidth=2) plot(xshelf,yshelf,'k--',linewidth=2) maketopo.theta_island = 260. (xisland,yisland) = latlong(1600.e3, maketopo.theta_island, 40., Rearth) x = linspace(xisland-1, xisland+1, 200) y = linspace(yisland-1, yisland+1, 200) X,Y = meshgrid(x,y) Z = maketopo.topo(X,Y) contour(X,Y,Z,[0,7,15]) axis('scaled') xlim([xisland-0.5, xisland+0.5]) ylim([yisland-0.5, yisland+0.5])