/[escript]/trunk/modellib/py_src/tsunami.py
ViewVC logotype

Contents of /trunk/modellib/py_src/tsunami.py

Parent Directory Parent Directory | Revision Log Revision Log

Revision 263 - (show annotations)
Tue Nov 29 09:55:16 2005 UTC (17 years, 3 months ago) by cochrane
File MIME type: text/x-python
File size: 47697 byte(s) 
Added visualisation of bathymmetry data, the coastline and the wave height
data as per earlier tsunami movies.  Also autogenerates the .params file for
ppmtompeg and generates the movie automatically.
1 #!/usr/bin/env python
2
3 # $Id$
4
5 import os, sys
6 import vtk
7 from esys.escript import *
8 from esys.escript.linearPDEs import LinearPDE
9 from esys.escript.modelframe import Model
10 import numarray
11 import math
12 import urllib
13
14 EPS=1.e-5
15
16 #====================================================================================================================================================
17 class GridData:
18 """
19 this object is used to store data on grid.
20 it will be replaced by Bruce at a later stage.
21
22 data[i,j] are data are x=j*s[0]+o[0] and y=i*s[1]+o[1]
23
24 for 0<=j<n[0] and 0<=i<n[1]
25 """
26 def __init__(self,s,o,n,data):
27 self.s=tuple(s)
28 self.o=tuple(o)
29 self.n=tuple(n)
30 self.data=data
31
32 def getOrigin(self):
33 return self.o
34 def getSpacing(self):
35 return self.s
36 def getData(self):
37 return self.data
38
39 def interpolate(self,x):
40 if hasattr(x,"convertToNumArray"):
41 x_array=x.convertToNumArray()
42 return_data_object=True
43 else:
44 x_array=numarray.array(x)
45 return_data_object=False
46 data=numarray.zeros(x_array.shape[0],numarray.Float)
47 ox,oy=self.getOrigin()
48 dx,dy=self.getSpacing()
49 data_array=self.getData()
50 i_dx=1
51 i_dy=1
52 for i in range(x_array.shape[0]):
53 x_long=x_array[i,0]-ox
54 x_lat=x_array[i,1]-oy
55 j0=min(max(int(x_long/dx),0),data_array.shape[1]-1-i_dy)
56 j1=min(max(int(x_lat/dy),0),data_array.shape[0]-1-i_dx)
57 f01=(x_long-j0*dx)/dx
58 f00=1.-f01
59 f11=(x_lat-j1*dy)/dy
60 f10=1.-f11
61 H00=data_array[j1,j0]
62 H01=data_array[j1,j0+i_dx]
63 H11=data_array[j1+i_dy,j0+i_dx]
64 H10=data_array[j1+i_dy,j0]
65 data[i]=(H00*f00+H01*f01)*f10+(H10*f00+H11*f01)*f11
66 if return_data_object:
67 out=Scalar(0,x.getFunctionSpace())
68 out.fillFromNumArray(data)
69 return out
70 else:
71 return data
72
73 def getVTK(self):
74 pass
75
76
77 class PointOnEarthSurface:
78 """
79 coordinates of a point on the surface of the Earth
80 """
81 def __init__(self,long=0,lat=0):
82 self.long=long
83 self.lat=lat
84 def __str__(self):
85 return "(%s,%s)"%(self.long,self.lat)
86
87 def __sub__(self,other):
88 return self.dist(other)
89
90 def split(self,p,t):
91 return PointOnEarthSurface(long=self.long+t*(p.long-self.long),lat=self.lat+t*(p.lat-self.lat))
92
93 def midPoint(self,other):
94 return PointOnEarthSurface(long=(self.long+other.long)/2,lat=(self.lat+other.lat)/2)
95
96 def dist(self,other):
97 return math.sqrt((self.long-other.long)**2+(self.lat-other.lat)**2)
98
99 class RegionOnEarthSurface:
100 """
101 defines an region by a south,east and north,west corner
102 """
103 RADIUS=6378.8e3
104 GRAVITY=9.81
105 WEST=0
106 NORTH=1
107 EAST=2
108 SOUTH=3
109 def __init__(self,west_south=PointOnEarthSurface(),east_north=PointOnEarthSurface(),resolution=1.):
110 if resolution<=0:
111 raise ValueError,"resolution must be positive"
112 if west_south.long>=east_north.long:
113 raise ValueError,"south west corner must be further east than west north corner"
114 if west_south.lat>=east_north.lat:
115 raise ValueError,"east south corner must be further down than west north corner"
116 if east_north.lat-west_south.lat < resolution/2:
117 raise ValueError,"latitude length of region must be 2*bigger then resolution"
118 if east_north.long-west_south.long < resolution/2:
119 raise ValueError,"longitude length of region must be 2*bigger then resolution"
120
121 self.west_south=west_south
122 self.east_north=east_north
123 self.resolution=resolution
124
125 def __str__(self):
126 return "RegionOnEarthSurface between %s and %s"%(str(self.west_south),str(self.east_north))
127
128 def isOnFace(self,p):
129 return self.isOnThisFace(p,self.SOUTH) or self.isOnThisFace(p,self.NORTH) or self.isOnThisFace(p,self.WEST) or self.isOnThisFace(p,self.EAST)
130 def isOnThisFace(self,p,face):
131 if face==self.WEST:
132 return self.west_south.long==p.long
133 if face==self.EAST:
134 return p.long==self.east_north.long
135 if face==self.SOUTH:
136 return self.west_south.lat==p.lat
137 if face==self.NORTH:
138 return p.lat==self.east_north.lat
139
140 def isBoxVertex(self,p):
141 return ( self.isOnThisFace(p,self.WEST) and self.isOnThisFace(p,self.SOUTH) ) or \
142 ( self.isOnThisFace(p,self.WEST) and self.isOnThisFace(p,self.NORTH) ) or \
143 ( self.isOnThisFace(p,self.EAST) and self.isOnThisFace(p,self.NORTH) ) or \
144 ( self.isOnThisFace(p,self.EAST) and self.isOnThisFace(p,self.SOUTH) )
145
146
147 def getFace(self,p):
148 # order is critical
149 if self.west_south.long==p.long: return self.WEST
150 if p.lat==self.east_north.lat: return self.NORTH
151 if p.long==self.east_north.long: return self.EAST
152 if self.west_south.lat==p.lat: return self.SOUTH
153
154 def comparePointsOnFace(self,p0,p1):
155 f0=self.getFace(p0)
156 f1=self.getFace(p1)
157
158 if f0<f1:
159 return -1
160 elif f0>f1:
161 return 1
162 else:
163 if f0 == self.WEST:
164 if p0.lat<p1.lat:
165 return -1
166 elif p0.lat>p1.lat:
167 return 1
168 else:
169 return 0
170 elif f0 == self.EAST:
171 if p0.lat<p1.lat:
172 return 1
173 elif p0.lat>p1.lat:
174 return -1
175 else:
176 return 0
177 elif f0 == self.NORTH:
178 if p0.long<p1.long:
179 return -1
180 elif p0.long>p1.long:
181 return 1
182 else:
183 return 0
184 else:
185 if p0.long<p1.long:
186 return 1
187 elif p0.long>p1.long:
188 return -1
189 else:
190 return 0
191
192
193 def isInRegion(self,p):
194 return self.west_south.long<=p.long \
195 and p.long<=self.east_north.long \
196 and self.west_south.lat<=p.lat \
197 and p.lat<=self.east_north.lat
198
199 def cutSegment(self,p0,p1):
200 t=None
201 p=None
202 tmp=self.interseptOfSegment(p0,p1,d=0,v=self.west_south.long)
203 if not tmp==None:
204 p_tmp=p0.split(p1,tmp)
205 if self.isInRegion(p_tmp) and t<tmp:
206 t=tmp
207 p=p_tmp
208
209 tmp=self.interseptOfSegment(p0,p1,d=0,v=self.east_north.long)
210 if not tmp==None:
211 p_tmp=p0.split(p1,tmp)
212 if self.isInRegion(p_tmp) and t<tmp:
213 t=tmp
214 p=p_tmp
215
216 tmp=self.interseptOfSegment(p0,p1,d=1,v=self.west_south.lat)
217 if not tmp==None:
218 p_tmp=p0.split(p1,tmp)
219 if self.isInRegion(p_tmp) and t<tmp:
220 t=tmp
221 p=p_tmp
222
223 tmp=self.interseptOfSegment(p0,p1,d=1,v=self.east_north.lat)
224 if not tmp==None:
225 p_tmp=p0.split(p1,tmp)
226 if self.isInRegion(p_tmp) and t<tmp:
227 t=tmp
228 p=p_tmp
229 return p
230
231 def interseptOfSegment(self,p0,p1,d=0,v=0.):
232 """
233 find t isn [0,1] such that p0+t*(p0-p1) cuts x[d]=v. if it does nit exist None is returned.
234 """
235 if d==0:
236 a=p0.long
237 b=p1.long
238 else:
239 a=p0.lat
240 b=p1.lat
241 if b==a:
242 if a==v:
243 t=0
244 else:
245 t=None
246 else:
247 t=(v-a)/(b-a)
248 if not (0<=t and t<=1): t=None
249 return t
250
251 class Polyline:
252 """
253 defines set of segments through a list of coordinates
254 """
255 def __init__(self,list_of_coordinates=[],name="none"):
256 c=[]
257 if len(list_of_coordinates)>0:
258 for i in range(len(list_of_coordinates)-1):
259 if list_of_coordinates[i]-list_of_coordinates[i+1]>0: c.append(list_of_coordinates[i])
260 c.append(list_of_coordinates[-1])
261 self.list_of_coordinates=c
262 self.name=name
263 def getDiameter(self):
264 out=0.
265 for p in self.list_of_coordinates:
266 for q in self.list_of_coordinates:
267 out=max(out,p-q)
268 return out
269 def isLoop(self):
270 if len(self)>0:
271 return not self[0]-self[-1]>EPS
272 else:
273 return True
274
275 def insert(self,index,coordinate):
276 """
277 insert point before index
278 """
279 if self.list_of_coordinates[index]-coordinate<EPS:
280 return index
281 elif self.list_of_coordinates[index-1]-coordinate<EPS:
282 return index-1
283 else:
284 self.list_of_coordinates.insert(index,coordinate)
285 return index
286
287 def split(self,index):
288 """
289 splits the polyline at point index
290 """
291 return Polyline(self.list_of_coordinates[:index],self.name),Polyline(self.list_of_coordinates[index:],self.name)
292
293 def __getitem__(self,s):
294 return self.list_of_coordinates.__getitem__(s)
295 def __iter__(self):
296 return iter(self.list_of_coordinates)
297
298 def __str__(self):
299 if self.isLoop():
300 out="loop["
301 else:
302 out="["
303 for i in self:
304 if out[-1]=="[":
305 out+="%s"%str(i)
306 else:
307 out+=",%s"%str(i)
308 return out+"]"
309
310 def __len__(self):
311 return len(self.list_of_coordinates)
312
313 def orientation(self):
314 """
315 returns the orientation of the polyline
316 """
317 if not self.isLoop():
318 raise TypeError,"polyline is not a loop"
319
320 integ=0.
321 for i in range(len(self.list_of_coordinates)-1):
322 p0=self.list_of_coordinates[i]
323 p1=self.list_of_coordinates[i+1]
324 integ+=(p1.lat-p0.lat)*(p1.long+p0.long)-(p1.long-p0.long)*(p1.lat+p0.lat)
325 if integ>=0.:
326 return 1.
327 else:
328 return -1.
329
330 def givePositiveOrientation(self):
331 if self.orientation()<0: self.list_of_coordinates.reverse()
332
333 class Coastline:
334 """
335 defines a coast line by a Polyline within a RegionOnEarthSurface
336 """
337 def __init__(self,region,name="none"):
338 self.region=region
339 self.name=name
340 self.polylines=[]
341 def __str__(self):
342 out=self.name+" in "+str(self.region)
343 for pl in self.polylines:
344 out+="\n"+str(pl)
345 return out
346 def makeTriangulation(self,west_south_is_water=True,east_south_is_water=True,west_north_is_water=True,east_north_is_water=True):
347 self.clean()
348 vertices=[]
349 segments=[]
350 holes=[]
351 vertices_on_face=[]
352 for pl in self.polylines:
353 if pl.getDiameter()>self.region.resolution:
354 short_pl=[pl[0]]
355 for i in range(1,len(pl)):
356 if pl[i]-short_pl[-1]>self.region.resolution/5:
357 short_pl.append(pl[i])
358 elif i==len(pl)-1:
359 short_pl[-1]=pl[i]
360 if pl.isLoop():
361 if len(short_pl)>3:
362 offset=len(vertices)
363 v_tmp=[short_pl[0]]
364 s_tmp=[]
365 for i in range(1,len(short_pl)):
366 if i==len(short_pl)-1:
367 s_tmp.append((len(v_tmp)-1+offset,offset))
368 else:
369 s_tmp.append((len(v_tmp)-1+offset,len(v_tmp)+offset))
370 v_tmp.append(short_pl[i])
371 vertices+=v_tmp
372 segments+=s_tmp
373 # find a point in the loop:
374 d_long=v_tmp[1].long-v_tmp[0].long
375 d_lat=v_tmp[1].lat-v_tmp[0].lat
376 l=math.sqrt(d_long**2+d_lat**2)
377 mid=v_tmp[0].midPoint(v_tmp[1])
378 n_long=-d_lat/l
379 n_lat=d_long/l
380 s=self.region.resolution
381 for seg in s_tmp:
382 p0=vertices[seg[0]]
383 p1=vertices[seg[1]]
384 a_long=p1.long-p0.long
385 a_lat=p1.lat-p0.lat
386 d=a_lat*n_long-a_long*n_lat
387 if abs(d)>0.:
388 t=((mid.lat-p0.lat)*n_long-(mid.long-p0.long)*n_lat)/d
389 if 0<=t and t<=1:
390 s_tmp=((mid.lat-p0.lat)*a_long-(mid.long-p0.long)*a_lat)/d
391 if s_tmp>EPS: s=min(s,s_tmp)
392 h=PointOnEarthSurface(long=mid.long+s/2*n_long,lat=mid.lat+s/2*n_lat)
393 holes.append(h)
394 else:
395 if len(short_pl)>1:
396 if self.region.isOnFace(short_pl[0]): vertices_on_face.append(short_pl[0])
397 if self.region.isOnFace(short_pl[-1]): vertices_on_face.append(short_pl[-1])
398 vertices.append(short_pl[0])
399 for i in range(1,len(short_pl)):
400 segments.append((len(vertices)-1,len(vertices)))
401 vertices.append(short_pl[i])
402 # put into the bounding box:
403 new_vertices=[]
404 if west_south_is_water:
405 new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long,lat=self.region.west_south.lat))
406 if east_south_is_water:
407 new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long,lat=self.region.west_south.lat))
408 if west_north_is_water:
409 new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long,lat=self.region.east_north.lat))
410 if east_north_is_water:
411 new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long,lat=self.region.east_north.lat))
412
413 # add new vertices if they don't exist yet
414 for q in new_vertices:
415 for p2 in vertices_on_face:
416 if p2-q<EPS:
417 q=None
418 raise ValueError,"coast line crosses boundary box vertex. This case is currenrly not supported."
419 if not q==None:
420 vertices.append(q)
421 vertices_on_face.append(q)
422 vertices_on_face.sort(self.region.comparePointsOnFace)
423 index=0
424 walking_on_water=west_south_is_water
425 l=len(vertices_on_face)
426 while index<l:
427 p1=vertices_on_face[(index+1)%l]
428 p0=vertices_on_face[index]
429 if walking_on_water:
430 segments.append((vertices.index(p0),vertices.index(p1)))
431 walking_on_water=False
432 else:
433 if self.region.isBoxVertex(p0):
434 segments.append((vertices.index(p0),vertices.index(p1)))
435 else:
436 walking_on_water=True
437 index+=1
438 return EarthTriangulation(vertices,segments,holes,self.region.resolution)
439
440 def clean(self):
441 """
442 cleans up the coast line by joining polylines to loops or connecting faces of the region
443 """
444 # find a poylines that are linked
445 while True:
446 k0=None
447 for pl in self.polylines:
448 if not pl.isLoop():
449 for k in [0,-1]:
450 for ql in self.polylines:
451 if not (ql==pl or ql.isLoop()):
452 for k2 in [0,-1]:
453 if ql[k2]-pl[k]<EPS:
454 pl0=pl
455 pl1=ql
456 k0=k
457 k1=k2
458 break
459 if not k0==None: break # ql
460 if not k0==None: break # k
461 if not k0==None: break # pl
462
463 if k0==None:
464 break
465 else:
466 self.polylines.remove(pl0)
467 self.polylines.remove(pl1)
468 pl0c=pl0.list_of_coordinates
469 pl1c=pl1.list_of_coordinates
470 if k0==0: pl0c.reverse()
471 if k1==-1: pl1c.reverse()
472 pl=Polyline(pl0c+pl1c[1:],pl0.name+" + "+pl1.name)
473 self.append(pl)
474
475 # find a polyline that is not a loop and has an end or start point not on the face of the region:
476 while True:
477 pl=None
478 k=None
479 for pl2 in self.polylines:
480 if not pl2.isLoop():
481 pl=pl2
482 if not self.region.isOnFace(pl[0]): k=0
483 if not self.region.isOnFace(pl[-1]): k=-1
484 if not k==None: break
485 if k==None: break
486 self.polylines.remove(pl)
487 d_min=50000.
488 k_min=None
489 for pl2 in self.polylines:
490 if not pl2.isLoop():
491 for k2 in [0,-1]:
492 if not self.region.isOnFace(pl2[k2]):
493 d2=pl2[k2]-pl[k]
494 if d2<d_min:
495 d_min=d2
496 pl_min=pl2
497 k_min=k2
498 if k_min==None:
499 raise ValueError,"cannot link coastline %s to any other coastline."%pl.name
500 plc=pl.list_of_coordinates
501 plc_min=pl_min.list_of_coordinates
502 if k==0: plc.reverse()
503 if k_min==-1: plc_min.reverse()
504 if d_min<EPS:
505 new_pl=Polyline(plc+plc_min[1:],pl.name+" + "+pl_min.name)
506 else:
507 new_pl=Polyline(plc+plc_min,pl.name+" + "+pl_min.name)
508 self.polylines.remove(pl_min)
509 self.append(new_pl)
510 # give positive orientation to loops:
511 for pl in self.polylines:
512 if pl.isLoop(): pl.givePositiveOrientation()
513
514 def append(self,polyline=Polyline()):
515 """append a polyline """
516 if len(polyline)>1:
517 pl=[]
518 outside_region=None
519 for i in range(len(polyline)):
520 if not self.region.isInRegion(polyline[i]):
521 outside_region=i
522 break
523 # pl.append(self.region.nudgeToFace(polyline[i]))
524 pl.append(polyline[i])
525 if not outside_region==None:
526 if outside_region==0:
527 for i in range(outside_region+1,len(polyline)):
528 if self.region.isInRegion(polyline[i]):
529 polyline.insert(i,self.region.cutSegment(polyline[i-1],\
530 polyline[i]))
531 pl1=polyline.split(i)[1]
532 self.append(pl1)
533 break
534 else:
535 # split polyline in two part first one is fully within the region the other starts with
536 # point outside the region
537 c=self.region.cutSegment(polyline[outside_region-1],polyline[outside_region])
538 i=polyline.insert(outside_region,c)
539 pl0,pl1=polyline.split(i+1)
540 self.append(pl0)
541 self.append(pl1)
542 else:
543 if len(pl)>1:
544 pply= Polyline(pl,polyline.name)
545 self.polylines.append(pply)
546
547 class EarthTriangulation:
548 GENERATOR="triangle -pqa%g %s"
549
550 def __init__(self,vertices=[],segments=[],holes=[],resolution=1.):
551 self.fn=os.tempnam()
552 # write triangle input file
553 poly_file=self.fn+".poly"
554 f=open(poly_file,"w")
555 f.writelines("%d %d %d %d\n"%(len(vertices),2,0,0))
556 for i in range(len(vertices)): f.writelines("%d %e %e\n"%(i,vertices[i].long,vertices[i].lat))
557 f.writelines("%d %d\n"%(len(segments),0))
558 for i in range(len(segments)): f.writelines("%d %d %d\n"%(i,segments[i][0],segments[i][1]))
559 f.writelines("%d\n"%(len(holes)))
560 for i in range(len(holes)): f.writelines("%d %e %e\n"%(i,holes[i].long,holes[i].lat))
561 f.close()
562 # start mesh generator:
563 os.system(self.GENERATOR%(resolution**2,poly_file))
564 # read mesh file:
565 self.node_coordinates=[]
566 self.node_tags=[]
567 self.node_ids=[]
568 self.triangles_nodes=[]
569 self.triangles_id=[]
570 node_file=open("%s.1.node"%self.fn,"r")
571 nodes=node_file.readline().strip().split()
572 nn=int(nodes[0])
573 for i in range(nn):
574 nodes=node_file.readline().strip().split()
575 self.node_coordinates.append((float(nodes[1]),float(nodes[2])))
576 self.node_tags.append(int(nodes[3]))
577 self.node_ids.append(int(nodes[0]))
578 node_file.close()
579 ele_file=open("%s.1.ele"%self.fn,"r")
580 elem=ele_file.readline().strip().split()
581 ne=int(elem[0])
582 for i in range(ne):
583 elem=ele_file.readline().strip().split()
584 self.triangles_id.append(int(elem[0]))
585 self.triangles_nodes.append((int(elem[1]),int(elem[2]),int(elem[3])))
586
587 ele_file.close()
588 # os.remove("%s.1.node"%self.fn)
589 # os.remove("%s.1.ele"%self.fn)
590 # os.remove("%s.ploy"%self.fn)
591
592 def getFinleyDomain(self):
593 from esys.finley import ReadMesh
594 finley_file=open("%s.msh"%self.fn,"w")
595 finley_file.writelines("%s\n2D-Nodes %d\n"%(self.fn,len(self.node_coordinates)))
596 for i in range(len(self.node_coordinates)):
597 finley_file.writelines("%s %s %s %e %e\n"%(self.node_ids[i],self.node_ids[i],self.node_tags[i],\
598 self.node_coordinates[i][0],self.node_coordinates[i][1]))
599
600 finley_file.writelines("Tri3 %d\n"%len(self.triangles_nodes))
601 for i in range(len(self.triangles_nodes)):
602 finley_file.writelines("%s 0 %s %s %s\n"%(self.triangles_id[i], \
603 self.triangles_nodes[i][0], \
604 self.triangles_nodes[i][1], \
605 self.triangles_nodes[i][2]))
606 finley_file.writelines("Line2 %d\n"%0)
607 finley_file.writelines("Line2_Contact %d\n"%0)
608 finley_file.writelines("Point1 %d\n"%0)
609 finley_file.close()
610 # get the mesh
611 out=ReadMesh("%s.msh"%self.fn)
612 # os.remove("%s.msh"%self.fn)
613 return out
614
615
616 #=================================
617 # Model interfaces:
618 #================================
619 class OceanRegionCollector(Model):
620 """
621
622
623 """
624 def __init__(self,debug=False):
625 Model.__init__(self,debug=debug)
626 self.declareParameter(coastline_source="http://jamboree.esscc.uq.edu.au/cgi-bin/doreen/datafile.txt?west=%%west%%&east=%%east%%&south=%%south%%&north=%%north%%&resolution=%%resolution%%&range360=%%range360%%&river=false&city=false&citytyp=wcity.dat&volcano=false&hotspot=false&shoreline=true&state=false&WSMdata=false&plate=false",
627 bathymetry_source="http://jamboree.esscc.uq.edu.au/cgi-bin/doreen/grdfile.xyz?west=%%west%%&east=%%east%%&south=%%south%%&north=%%north%%&resolution=%%resolution%%",
628 resolution=1.,
629 south=0.,
630 north=10.,
631 east=0.,
632 west=20.,
633 range360=True,
634 coastline_stream=None,
635 bathymetry_stream=None)
636
637
638 def doInitialization(self):
639 """
640 Initializes the ocean region
641 """
642 c=self.__mergeParameters(self.coastline_source)
643 b=self.__mergeParameters(self.bathymetry_source)
644 self.coastline_stream=urllib.urlopen(c)
645 self.bathymetry_stream=urllib.urlopen(b)
646
647 def __mergeParameters(self,txt):
648 return txt.replace("%%west%%",str(self.west))\
649 .replace("%%east%%",str(self.east))\
650 .replace("%%south%%",str(self.south))\
651 .replace("%%north%%",str(self.north))\
652 .replace("%%resolution%%",str(self.resolution)) \
653 .replace("%%range360%%",str(self.range360).lower())
654
655 class Bathymetry(Model):
656 """
657 generates the bathymetry data within a region on the earth
658 """
659 def __init__(self,debug=False):
660 Model.__init__(self,debug=debug)
661 self.declareParameter(source="none",
662 bathymetry=1.)
663
664 def doInitialization(self):
665 """
666 Initializes the
667 """
668 if hasattr(self.source,"readline"):
669 f=self.source
670 else:
671 f=open(filename,"r")
672 x_grd_list=[]
673 y_grd_list=[]
674 data_grd_list=[]
675 line=f.readline().strip()
676 while line!="":
677 v=line.split()
678 x_grd_list.append(float(v[0]))
679 y_grd_list.append(float(v[1]))
680 data_grd_list.append(float(v[2]))
681 line=f.readline().strip()
682 self.trace("%s data have been read from %s."%(len(data_grd_list),self.source))
683 data_grd=numarray.array(data_grd_list)
684 x_grd=numarray.array(x_grd_list)
685 y_grd=numarray.array(y_grd_list)
686 if len(x_grd)<2:
687 raise ValueError,"%s: data base is two small"%str(self)
688 ox=x_grd[0]
689 oy=y_grd[0]
690 diam=max(abs(x_grd[len(x_grd)-1]-ox),abs(y_grd[len(y_grd)-1]-oy))
691 dx=x_grd[1]-ox
692 nx=1
693 nx=1
694 while abs(x_grd[nx]-ox)>1.e-10*diam:
695 nx+=1
696 dy=y_grd[nx]-oy
697 ny=len(x_grd)/nx
698 data_grd.resize((ny,nx))
699 self.bathymetry=GridData(s=[dx,dy],o=[ox,oy],n=[nx,ny],data=data_grd)
700 self.trace("%s x %s grid with %s x %s spacing."%(nx,ny,dx,dy))
701
702
703 class OceanRegion(Model):
704 """
705 generates the ocean region with a coast line and a bathymetry
706
707 """
708 def __init__(self,debug=False):
709 Model.__init__(self,debug=debug)
710 self.declareParameter(domain=None, \
711 resolution=1.,
712 south=0.,
713 north=10.,
714 east=0.,
715 west=20.,
716 bathymetry=None,
717 bathymetry_data=None,
718 coastline=None,
719 source="none")
720
721 def doInitialization(self):
722 """
723 Initializes the ocean region
724 """
725 if hasattr(self.source,"readline"):
726 f=self.source
727 data_name=f.geturl()
728 else:
729 f=open(self.source,"r")
730 data_name=self.source
731
732 segs=[]
733 name=""
734 line=f.readline()
735 while not line=="":
736 line=line.strip()
737 if line[:7]=="#region":
738 data=line[line.index("[")+1:line.index("]")].split(",")
739 reg=RegionOnEarthSurface(PointOnEarthSurface(lat=self.south,long=self.west),PointOnEarthSurface(lat=self.north,long=self.east),self.resolution)
740 self.trace(str(reg))
741 self.coastline=Coastline(region=reg,name=data_name)
742 elif line.find("Shore Bin")>-1:
743 self.coastline.append(Polyline(segs,name))
744 name=line[2:]
745 segs=[]
746 if not (line=="" or line[0]=="#" or line[0]==">") :
747 x=line.split()
748 segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
749 line=f.readline()
750 self.coastline.append(Polyline(segs,name))
751 d=self.bathymetry_data.interpolate([[self.east,self.south],[self.west,self.south],[self.east,self.north],[self.west,self.north]])
752 self.domain=self.coastline.makeTriangulation(east_south_is_water=d[0]<=0,west_south_is_water=d[1]<=0,east_north_is_water=d[2]<=0,west_north_is_water=d[3]<=0).getFinleyDomain()
753 self.bathymetry=maximum(-self.bathymetry_data.interpolate(Function(self.domain).getX()),0.)
754
755
756 class TsunamiSource(Model):
757 """
758 defines a wave in Gaussean form between start and end.
759 """
760 GAMMA=0.05
761 def __init__(self,debug=False):
762 Model.__init__(self,debug=debug)
763 self.declareParameter(domain=None,
764 start_lat=-10.,
765 start_long=110.,
766 end_lat=-12.,
767 end_long=120.,
768 width=5.,
769 decay_zone=0.1,
770 amplitude=1.,
771 wave_height=1.)
772
773 def doInitialization(self):
774 """
775 set initial wave form
776 """
777 beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone
778 x=self.domain.getX()
779 x_long=x[0]
780 x_lat=x[1]
781 mid_long=(self.start_long+self.end_long)/2
782 mid_lat=(self.start_lat+self.end_lat)/2
783 dist=math.sqrt((mid_long-self.end_long)**2+(mid_lat-self.end_lat)**2)
784 a=(self.start_long-mid_long)/dist
785 b=(self.start_lat-mid_lat)/dist
786 self.trace("source length = %s"%(dist*2.))
787 x_long_hat= a*(x_long-mid_long)+b*(x_lat-mid_lat)
788 x_lat_hat=-b*(x_long-mid_long)+a*(x_lat-mid_lat)
789 # x_lat-direction
790 s=abs(x_lat_hat)-self.width
791 m=s.whereNegative()
792 f1=(1.-m)*exp(-(s*beta)**2)+m
793 # x_long-direction
794 s=abs(x_long_hat)-dist
795 m=s.whereNegative()
796 f0=(1.-m)*exp(-(s*beta)**2)+m
797 self.wave_height=f1*f0*self.amplitude
798
799 #====================================================================================================================================================
800 class TsunamiInDeepWater(Model):
801 """
802 Runs the deep water tsunami model based on a simplfied form of the shallow water equation.
803
804
805 M{d^2 h/dt^2 =div(c grad(h)) }
806
807 where h is the wave height above sea level, and with H the depth of the water level and g is the gravity constant
808 c=sqrt(g*H).
809
810 The simulation uses the Verlet scheme.
811
812 """
813 def __init__(self,debug=False):
814 Model.__init__(self,debug=debug)
815 self.declareParameter(domain=None, \
816 wave_height=1.,
817 wave_velocity=0.,
818 initial_time_step=None,
819 bathymetry=1.,
820 safety_factor=1.)
821
822 def doInitialization(self):
823 """
824 Initializes the time integartion scheme
825 """
826 self.__pde=LinearPDE(self.domain)
827 self.__pde.setSolverMethod(self.__pde.LUMPING)
828 self.__pde.setValue(D=1.)
829 self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/RegionOnEarthSurface.RADIUS**2
830 c_max=math.sqrt(Lsup(self.__c2))
831 self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))
832 if self.initial_time_step==None: self.initial_time_step=self.__dt
833 self.trace("maximum wave velocity %s m/sec"%c_max)
834 self.trace("Time step size is %s sec"%self.__dt)
835
836
837 def getSafeTimeStepSize(self,dt):
838 """
839 returns new step size
840
841 @param dt: last time step size used
842 @type dt: C{float}
843 @return: time step size that can savely be used
844 @rtype: C{float}
845 """
846 return self.__dt
847
848 def doStepPostprocessing(self,dt):
849 """
850 perform the time step using the Valet scheme
851
852 @param dt: time step size to be used
853 @type dt: C{float}
854 """
855 self.__pde.setValue(X=-self.__c2*grad(self.wave_height))
856
857 new_height=self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()
858
859 self.wave_velocity=(new_height-self.wave_height)/dt
860 self.wave_height=new_height
861 self.initial_time_step=dt
862 self.trace("Wave height range is %e %e"%(inf(self.wave_height),sup(self.wave_height)))
863
864 class SurfMovie(Model):
865 """
866 movie from a wave propagation on the sea
867
868 @ivar time: current time
869 @ivar bathymetry: scalar data set
870 @ivar wave_height: vector data set
871 @ivar filename: name of the movie file
872 """
873 def __init__(self,debug=False):
874 Model.__init__(self,debug=debug)
875 self.declareParameter(bathymetry=1.,
876 wave_height=1.,
877 coastline=None,
878 t=0.,
879 dt=1.,
880 south=2.,
881 north=5.,
882 east=3.,
883 west=15.,
884 filename="movie.mpg")
885
886 self.paramsFileString = "REFERENCE_FRAME DECODED\n"
887 self.paramsFileString += "FRAME_RATE 24\n"
888 self.paramsFileString += "OUTPUT %s\n" % self.filename
889 self.paramsFileString += "PATTERN IBBPBBPBB\n"
890 self.paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
891 self.paramsFileString += "GOP_SIZE 20\n"
892 self.paramsFileString += "BSEARCH_ALG CROSS2\n"
893 self.paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
894 self.paramsFileString += "IQSCALE 10\n"
895 self.paramsFileString += "PQSCALE 11\n"
896 self.paramsFileString += "BQSCALE 16\n"
897 self.paramsFileString += "RANGE 8\n"
898 self.paramsFileString += "SLICES_PER_FRAME 1\n"
899 self.paramsFileString += "BASE_FILE_FORMAT PNM\n"
900 self.paramsFileString += "INPUT_DIR \n"
901 self.paramsFileString += "INPUT_CONVERT *\n"
902 self.paramsFileString += "INPUT\n"
903
904 self.firstFrame = True
905
906 def doInitialization(self):
907 """
908 Initializes the time integration scheme
909 """
910 self.__fn=os.tempnam()+".xml"
911 self.__frame_name=os.tempnam()
912 self.__next_t=self.dt
913 # self.coastline.getVTK()
914 # self.bathymetry.getVTK()
915 # wndow(south,west,north,east)
916
917 # the bathymmetry colourmap
918 data = []
919 data.append([-8000, 0, 0, 0])
920 data.append([-7000, 0, 5, 25])
921 data.append([-6000, 0, 10, 50])
922 data.append([-5000, 0, 80, 125])
923 data.append([-4000, 0, 150, 200])
924 data.append([-3000, 86, 197, 184])
925 data.append([-2000, 172, 245, 168])
926 data.append([-1000, 211, 250, 211])
927 data.append([0, 16, 48, 16])
928 data.append([300, 70, 150, 50])
929 data.append([500, 146, 126, 60])
930 data.append([1000, 198, 178, 80])
931 data.append([1250, 250, 230, 100])
932 data.append([1500, 250, 234, 126])
933 data.append([1750, 252, 238, 152])
934 data.append([2000, 252, 243, 177])
935 data.append([2250, 253, 249, 216])
936 data.append([2500, 255, 255, 255])
937
938 # the amount to scale the data by
939 scale = 255.0
940 numColours = len(data)
941
942 # convert the colourmap into something vtk is more happy with
943 height = numarray.zeros(numColours, numarray.Float)
944 red = numarray.zeros(numColours, numarray.Float)
945 green = numarray.zeros(numColours, numarray.Float)
946 blue = numarray.zeros(numColours, numarray.Float)
947 for i in range(numColours):
948 (h, r, g, b) = data[i]
949 height[i] = float(h)
950 red[i] = float(r)/scale
951 green[i] = float(g)/scale
952 blue[i] = float(b)/scale
953
954 print "Loading bathymmetry data..."
955 # grab the z data
956 bathZData = self.bathymetry.getData()
957
958 # get the origin
959 bathOrigin = self.bathymetry.getOrigin()
960 # get the delta_x and delta_y
961 bathSpacing = self.bathymetry.getSpacing()
962
963 # now construct the x and y data from the spacing and the origin
964 numXPoints = bathZData.shape[1]
965 numYPoints = bathZData.shape[0]
966 numPoints = numXPoints*numYPoints
967
968 bathXData = numarray.zeros(numXPoints, numarray.Float)
969 bathYData = numarray.zeros(numYPoints, numarray.Float)
970 for i in range(numXPoints):
971 bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
972
973 for i in range(numYPoints):
974 bathYData[i] = bathOrigin[1] + i*bathSpacing[1]
975
976 # calculate the appropriate window size
977 dataWidth = max(bathXData) - min(bathXData)
978 dataHeight = max(bathYData) - min(bathYData)
979 winHeight = 600
980 winWidth = int(dataWidth*float(winHeight)/dataHeight)
981
982 print "Done loading bathymmetry data"
983
984 ### now do the vtk stuff
985
986 # make sure rendering is offscreen
987 factGraphics = vtk.vtkGraphicsFactory()
988 factGraphics.SetUseMesaClasses(1)
989
990 factImage = vtk.vtkImagingFactory()
991 factImage.SetUseMesaClasses(1)
992
993 # make the bathymmetry colourmap
994 transFunc = vtk.vtkColorTransferFunction()
995 transFunc.SetColorSpaceToRGB()
996 for i in range(numColours):
997 transFunc.AddRGBPoint(height[i], red[i], green[i], blue[i])
998 h = height[i]
999 while i < numColours-1 and h < height[i+1]:
1000 h += 1
1001 transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1002
1003 # set up the lookup table for the wave data
1004 refLut = vtk.vtkLookupTable()
1005 self.lutTrans = vtk.vtkLookupTable()
1006 refLut.Build()
1007 alpha = 0.4 # alpha channel value
1008 for i in range(256):
1009 (r,g,b,a) = refLut.GetTableValue(255-i)
1010 if g == 1.0 and b < 0.5 and r < 0.5:
1011 self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1012 else:
1013 self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1014
1015 print "Generating the bathymmetry vtk object..."
1016
1017 # set up the renderer and the render window
1018 self.ren = vtk.vtkRenderer()
1019 self.renWin = vtk.vtkRenderWindow()
1020 self.renWin.AddRenderer(self.ren)
1021 self.renWin.SetSize(winWidth, winHeight)
1022 self.renWin.OffScreenRenderingOn()
1023
1024 # make an unstructured grid
1025 bathGrid = vtk.vtkUnstructuredGrid()
1026
1027 # make the points for the vtk data
1028 bathPoints = vtk.vtkPoints()
1029 bathPoints.SetNumberOfPoints(numPoints)
1030
1031 # make the vtk bathymmetry data set
1032 bathData = vtk.vtkFloatArray()
1033 bathData.SetNumberOfComponents(1)
1034 bathData.SetNumberOfTuples(numPoints)
1035
1036 # add the points and data
1037 count = 0
1038 for i in range(numXPoints):
1039 for j in range(numYPoints):
1040 bathPoints.InsertPoint(count, bathXData[i], bathYData[j], 0.0)
1041 bathData.InsertTuple1(count, bathZData[j,i])
1042 count += 1
1043
1044 # set the data to the grid
1045 bathGrid.SetPoints(bathPoints)
1046 bathGrid.GetPointData().SetScalars(bathData)
1047
1048 # do a delaunay on the grid
1049 bathDelaunay = vtk.vtkDelaunay2D()
1050 bathDelaunay.SetInput(bathGrid)
1051 bathDelaunay.SetTolerance(0.001)
1052 bathDelaunay.SetAlpha(2.5)
1053
1054 # set up the mapper
1055 bathMapper = vtk.vtkDataSetMapper()
1056 bathMapper.SetInput(bathDelaunay.GetOutput())
1057 bathMapper.SetLookupTable(transFunc)
1058
1059 # set up the actor
1060 bathActor = vtk.vtkActor()
1061 bathActor.SetMapper(bathMapper)
1062
1063 self.ren.AddActor(bathActor)
1064
1065 ### now add the coastline
1066 print "Loading the coastline data..."
1067
1068 # make the coastline points
1069 coastPoints = vtk.vtkPoints()
1070
1071 # make the coastline grid
1072 coastGrid = vtk.vtkUnstructuredGrid()
1073
1074 # now point the points and lines into the grid
1075 totalCoastPoints = 0
1076 for polyline in self.coastline.polylines:
1077 numPoints = len(polyline)
1078 coastLine = vtk.vtkPolyLine()
1079 coastLine.GetPointIds().SetNumberOfIds(numPoints)
1080 j = 0
1081 for point in polyline:
1082 coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1083 coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1084 j += 1
1085 coastGrid.InsertNextCell(coastLine.GetCellType(),
1086 coastLine.GetPointIds())
1087 totalCoastPoints += numPoints
1088
1089 coastGrid.SetPoints(coastPoints)
1090
1091 # make the coast's mapper
1092 coastMapper = vtk.vtkDataSetMapper()
1093 coastMapper.SetInput(coastGrid)
1094
1095 # make its actor
1096 coastActor = vtk.vtkActor()
1097 coastActor.SetMapper(coastMapper)
1098 coastActor.GetProperty().SetColor(0,0,0)
1099
1100 # add the actor to the renderer
1101 self.ren.AddActor(coastActor)
1102
1103 # set up the actor for the wave
1104 self.waveActor = vtk.vtkActor()
1105
1106 # add the actor to the renderer
1107 self.ren.AddActor(self.waveActor)
1108
1109 print "Done loading the coastline data"
1110
1111 def doStepPostprocessing(self, dt):
1112 """
1113 Does any necessary postprocessing after each step
1114
1115 @param dt:
1116 """
1117 if self.t>=self.__next_t:
1118 print self.t,"write ",Lsup(self.wave_height)
1119 saveVTK(self.__fn,h=self.wave_height)
1120 # vtkobj=...
1121 # save(self.__frame_name)
1122
1123 # get the wave data
1124 waveDomain = self.wave_height.getDomain().getX()
1125 waveX = waveDomain[0].convertToNumArray()
1126 waveY = waveDomain[1].convertToNumArray()
1127 waveZ = self.wave_height.convertToNumArray()
1128
1129 numPoints = len(waveZ)
1130
1131 # make the points
1132 wavePoints = vtk.vtkPoints()
1133 wavePoints.SetNumberOfPoints(numPoints)
1134
1135 # make the vtk data array
1136 waveData = vtk.vtkFloatArray()
1137 waveData.SetNumberOfComponents(1)
1138 waveData.SetNumberOfTuples(numPoints)
1139 waveData.SetName("data")
1140
1141 # put the data into the points and array
1142 for i in range(numPoints):
1143 wavePoints.InsertPoint(i, waveX[i], waveY[i], 0.0)
1144 waveData.InsertTuple1(i, waveZ[i])
1145
1146 # make the grid
1147 waveGrid = vtk.vtkUnstructuredGrid()
1148 waveGrid.SetPoints(wavePoints)
1149 waveGrid.GetPointData().AddArray(waveData)
1150 waveGrid.GetPointData().SetActiveScalars("data")
1151
1152 (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1153 print "Wave height range %f - %f" % (zMin, zMax)
1154
1155 # do a delaunay on the data
1156 waveDelaunay = vtk.vtkDelaunay2D()
1157 waveDelaunay.SetInput(waveGrid)
1158 waveDelaunay.SetTolerance(0.001)
1159 waveDelaunay.SetAlpha(2.5)
1160
1161 # make a mapper for the grid
1162 waveMapper = vtk.vtkDataSetMapper()
1163 waveMapper.SetInput(waveDelaunay.GetOutput())
1164 waveMapper.SetLookupTable(self.lutTrans)
1165 waveMapper.SetScalarRange(zMin, zMax)
1166
1167 self.waveActor.SetMapper(waveMapper)
1168
1169 # do stuff here that can only be done on the first frame
1170 if self.firstFrame:
1171 # zoom in a bit
1172 self.ren.GetActiveCamera().Zoom(2.0)
1173 self.firstFrame = False
1174
1175 # render the window
1176 self.renWin.Render()
1177
1178 # convert the render window to an image
1179 win2imgFilter = vtk.vtkWindowToImageFilter()
1180 win2imgFilter.SetInput(self.renWin)
1181
1182 # the image name
1183 imgFname = "%s%f.pnm" % (self.__frame_name, self.__next_t)
1184
1185 # save the image to file
1186 outWriter = vtk.vtkPNMWriter()
1187 outWriter.SetInput(win2imgFilter.GetOutput())
1188 outWriter.SetFileName(imgFname)
1189 outWriter.Write()
1190 print "Wrote %s" % imgFname
1191 os.system("display %s" % imgFname)
1192 self.paramsFileString += "%s\n" % imgFname
1193
1194 self.__next_t+=self.dt
1195
1196 def getSafeTimeStepSize(self,dt):
1197 """
1198 returns new step size
1199
1200 @param dt: last time step size used
1201 @type dt: C{float}
1202 @return: time step size that can savely be used
1203 @rtype: C{float}
1204 """
1205 return self.__next_t-self.t
1206
1207 def doFinalization(self):
1208 """
1209 Finalises the visualisation. For instance, makes a movie of the
1210 image files.
1211 """
1212 # make the movie into self.filename
1213 self.paramsFileString += "END_INPUT\n"
1214 self.paramsFileString += "PIXEL HALF\n"
1215 self.paramsFileString += "ASPECT_RATIO 1\n"
1216
1217 # write the params string to file
1218 fp = open("%s.params" % self.filename, "w")
1219 fp.write(self.paramsFileString + '\n')
1220 fp.close()
1221
1222 print "Performing conversion to mpeg"
1223 convertString = "ppmtompeg %s.params" % self.filename
1224 result = os.system(convertString)
1225 if result != 0:
1226 print "An error occurred in mpeg conversion"
1227
1228 if __name__=="__main__":
1229 from esys.escript.modelframe import Link,Simulation
1230 from esys.modellib.input import Sequencer
1231
1232 oc=OceanRegionCollector()
1233 oc.resolution=2
1234 oc.south=-45.5
1235 oc.north=-5.4
1236 oc.east=161.1
1237 oc.west=108.2
1238 oc.range360=True
1239
1240
1241 b=Bathymetry()
1242 b.source=Link(oc,"bathymetry_stream")
1243
1244 oreg=OceanRegion()
1245 oreg.source=Link(oc,"coastline_stream")
1246 oreg.resolution=Link(oc,"resolution")
1247 oreg.south=Link(oc,"south")
1248 oreg.north=Link(oc,"north")
1249 oreg.east=Link(oc,"east")
1250 oreg.west=Link(oc,"west")
1251 oreg.bathymetry_data=Link(b,"bathymetry")
1252
1253 src=TsunamiSource()
1254 src.domain=Link(oreg,"domain")
1255 src.start_lat=-10.
1256 src.end_lat=-12.
1257 src.start_long=110.
1258 src.end_long=120.
1259 src.width=0.1
1260 src.decay_zone=0.01
1261 src.amplitude=1.
1262
1263 ts=TsunamiInDeepWater()
1264 ts.domain=Link(oreg,"domain")
1265 ts.wave_height=Link(src,"wave_height")
1266 ts.wave_velocity=0.
1267 ts.bathymetry=Link(oreg,"bathymetry")
1268
1269 sq=Sequencer()
1270 sq.t_end=100000.
1271
1272 sm=SurfMovie()
1273 sm.bathymetry=Link(b,"bathymetry")
1274 sm.wave_height=Link(ts,"wave_height")
1275 sm.coastline=Link(oreg,"coastline")
1276 sm.t=Link(sq,"t")
1277 sm.dt=5000.
1278 sm.filename="movie.mpg"
1279
1280 s=Simulation([sq,oc,b,oreg,src,ts,sm])
1281 # s.writeXML()
1282 s.run()
1283
1284 # vim: expandtab shiftwidth=4:
  ViewVC Help
Powered by ViewVC 1.1.26