Changeset 5419

Timestamp:

Jun 24, 2008, 12:37:05 PM (17 years ago)

Author:

jakeman

Message:

adjustedproject file to read in order of boundary nodes from file

File:

• anuga_work/development/boxingday08/project.py (modified) (6 diffs)

anuga_work/development/boxingday08/project.py

@@ -2 +2 @@
 from Scientific.IO.NetCDF import NetCDFFile
 from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
-     compress,zeros,fabs
+    compress,zeros,fabs,allclose
 from anuga.utilities.polygon import inside_polygon,read_polygon
 from os import sep
 from time import localtime, strftime, gmtime
-def create_sts_boundary(filename,verbose=False):
-    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
+
+def create_sts_boundary_order_file(order_filename,stsfilename,verbose=False):
+    """
+    Note This is not a robust algorithm. Be Careful when applying that 
+    polygon does not intersect itself and that resulting boundary actually
+    follows the path it is supposed too.
+    """
+    fid = NetCDFFile(stsfilename+'.sts', 'r')    #Open existing file for read
     x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
     y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
@@ -33 +39 @@
     boundary[0][1]=y[west_index]
     indices.remove(west_index)
+    order=[]
     for i in range(n-1):
          min_dist=1e12
@@ -46 +53 @@
          boundary[i+1][0]=x[index]
          boundary[i+1][1]=y[index]
+         order.append(index)
     #ensure boundary is a polygon i.e no intesections
     #e.g.
@@ -52 +60 @@
               
     #So far i cannot tell if the wrong point is closer.
+    d=","
+    fid=open(order_filename,'w')
+    header="index,longitude,latitude\n"
+    fid.write(header)
+    line=str(west_index)+d+str(x[west_index])+d+\
+        str(y[west_index])+"\n"
+    fid.write(line)
+    for i in order:
+            line=str(i)+d+str(x[i])+d+\
+                str(y[i])+"\n"
+            fid.write(line)
+    fid.close()
     return boundary
+            
 
 #read in polyline boundary
-from anuga.shallow_water.data_manager import urs2sts
+from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
 base_name='polyline'
 import os
 if os.path.exists(base_name+'.sts'):
-        print 'sts boundary file already created.'
+        print 'sts boundary file already created.'
 else:
-    urs2sts(base_name,mean_stage=0.0,verbose=False)
-boundary=create_sts_boundary(base_name,verbose=True)
-urs_boundary=boundary.tolist()
+        urs2sts(base_name,mean_stage=0.0,verbose=False)
+
+order_filename='sts_order.txt'
+boundary=create_sts_boundary_order_file(order_filename,base_name,verbose=False)
+urs_boundary=create_sts_boundary(order_filename,base_name,lat_long=False)
+assert allclose(boundary.tolist(),urs_boundary)
 
 #Read in large domain boundary
@@ -70 +94 @@
 bp = read_polygon(dir+sep+extent)
 
+#Now order of boundary is correct clip sts boundary to small region for
+#this particular scenario
 indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False)
 bounding_polygon=[]
@@ -76 +102 @@
     bounding_polygon.append(urs_boundary[i])
 
+#Append clipped sts boundary to large domain boundary
 bounding_polygon.extend(bp[2:len(bp)])