Changeset 1661

Timestamp:

Aug 1, 2005, 12:11:53 PM (19 years ago)

Author:

ole

Message:

Brought xya2pts back

File:

• inundation/ga/storm_surge/pyvolution/data_manager.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -51 +51 @@
 from Numeric import concatenate
 
-from coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE
+from coordinate_transforms.geo_reference import Geo_reference
 
 def make_filename(s):
@@ -921 +921 @@
     return cls(domain, mode)
 
+def xya2pts(basename_in, basename_out=None, verbose=False,
+            #easting_min=None, easting_max=None,
+            #northing_min=None, northing_max=None,
+            attribute_name = 'elevation'):
+    """Read points data from ascii (.xya)
+
+    Example:
+
+              x(m)        y(m)        z(m)
+         0.00000e+00  0.00000e+00  1.3535000e-01
+         0.00000e+00  1.40000e-02  1.3535000e-01
+
+    
+
+    Convert to NetCDF pts format which is
+
+    points:  (Nx2) Float array
+    elevation: N Float array
+
+    Only lines that contain three numeric values are processed
+    """
+
+    import os 
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import Float, arrayrange, concatenate
+
+    root, ext = os.path.splitext(basename_in)
+
+    if ext == '': ext = '.xya'
+
+    #Get NetCDF
+    infile = open(root + ext, 'r')  #Open existing xya file for read
+
+    if verbose: print 'Reading xya points from %s' %(root + ext)
+
+    points = []
+    attribute = []
+    for i, line in enumerate(infile.readlines()):
+        fields = line.split()
+
+        try:
+            assert len(fields) == 3
+        except:    
+            print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line) 
+
+        try:
+            x = float( fields[0] )
+            y = float( fields[1] )
+            z = float( fields[2] )            
+        except:
+            continue
+
+        points.append( [x, y] )
+        attribute.append(z)
+
+
+    #Get output file
+    if basename_out == None:
+        ptsname = root + '.pts'
+    else:
+        ptsname = basename_out + '.pts'
+
+    if verbose: print 'Store to NetCDF file %s' %ptsname
+    # NetCDF file definition
+    outfile = NetCDFFile(ptsname, 'w')
+
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'NetCDF pts format for compact and '\
+                          'portable storage of spatial point data'
+    
+
+    #Georeferencing
+    from coordinate_transforms.geo_reference import Geo_reference
+    Geo_reference().write_NetCDF(outfile)
+    
+
+    outfile.createDimension('number_of_points', len(points))
+    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+
+    # variable definitions
+    outfile.createVariable('points', Float, ('number_of_points',
+                                             'number_of_dimensions'))
+    outfile.createVariable(attribute_name, Float, ('number_of_points',))
+
+    # Get handles to the variables
+    nc_points = outfile.variables['points']
+    nc_attribute = outfile.variables[attribute_name]
+
+    #Store data
+    nc_points[:, :] = points
+    nc_attribute[:] = attribute
+
+    infile.close()
+    outfile.close()
+    
 def dem2pts(basename_in, basename_out=None, verbose=False,
             easting_min=None, easting_max=None,
@@ -1144 +1240 @@
     if origin is None:
 
-        # get geo_reference
+        #Get geo_reference
         #sww files don't have to have a geo_ref
         try:
             geo_reference = Geo_reference(NetCDFObject=fid)
         except AttributeError, e:
-            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
+            geo_reference = Geo_reference() #Default georef object
 
         xllcorner = geo_reference.get_xllcorner()