Changeset 5581 for anuga_work/production/perth/get_gauges.py

Timestamp:

Jul 28, 2008, 6:20:17 PM (16 years ago)

Author:

sexton

Message:

fix script to get gauges from sts file and test run_perth with tide=0

File:

• anuga_work/production/perth/get_gauges.py (modified) (4 diffs)

anuga_work/production/perth/get_gauges.py

@@ -9 +9 @@
 from anuga.utilities.numerical_tools import ensure_numeric
 from Scientific.IO.NetCDF import NetCDFFile
-
 
 def get_sts_gauge_data(filename,polygon,verbose=False):
@@ -31 +30 @@
         #show()
 
+    
     basename='sts_gauge'
     quantity_names=['stage','xmomentum','ymomentum']
@@ -38 +38 @@
         if inside_points is not None:
             quantities[name] = take(quantities[name],indices,1)
-            for j in range(len(x)-1):
-                locx=int(x[j])
-                locy=int(y[j])
-                fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
-                s = 'time, stage, xmomentum, ymomentum \n'
-                fid_sts.write(s)
-                stage = quantities[name][:,j]
-                xmomentum = quantities[name][:,j]
-                ymomentum = quantities[name][:,j]
-                for k in range(len(time)-1):
-                    s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
-                    fid_sts.write(s)
         else:
             msg = 'No gauges found inside polygon'
             raise msg
-            fid_sts.close()
+            
+    for j in range(len(x)-1):
+        locx=int(x[j])
+        locy=int(y[j])
+        fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
+        s = 'time, stage, xmomentum, ymomentum \n'
+        fid_sts.write(s)
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+        for k in range(len(time)-1):
+            s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
+            fid_sts.write(s)
+
+        fid_sts.close()      
+
     fid.close()
-
+        
     return inside_points,quantities,elevation,time
 
@@ -61 +64 @@
 points,quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),polygon,verbose=False)
 
-
 assert len(points)==len(elevation)==len(quantities['stage'][0,:])