Changeset 928 for inundation/ga/storm_surge/pyvolution/data_manager.py

Timestamp:

Feb 18, 2005, 5:15:46 PM (20 years ago)

Author:

ole

Message:

Incorporated variable false eastings and northings (just in case:-).

File:

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

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -872 +872 @@
     dem_elevation = infile.variables['elevation']
 
+    zone = infile.zone[0]
+    false_easting = infile.false_easting[0]
+    false_northing = infile.false_northing[0]
+
+    #Text strings
+    projection = infile.projection
+    datum = infile.datum
+    units = infile.units
+    
+
     #Get output file                                
     xyaname = root + '.pts'
@@ -884 +894 @@
 
     #Georeferencing
-    outfile.zone = 56   #FIXME: Must be read from somewhere. Talk to Don.
+    outfile.zone = zone
     outfile.xllcorner = xllcorner #Easting of lower left corner
     outfile.yllcorner = yllcorner #Northing of lower left corner
-    
+    outfile.false_easting = false_easting
+    outfile.false_northing =false_northing
+    
+    outfile.projection = projection
+    outfile.datum = datum
+    outfile.units = units
+
 
     #Grid info (FIXME: probably not going to be used, but heck)
@@ -938 +954 @@
 
     Convert to NetCDF format (.dem) mimcing the ASCII format closely.
+
+
+    An accompanying file with same basename but extension .prj must exist
+    and is used to fix the UTM zone, datum, false northings and eastings.
+
+    The prj format is assumed to be as
+    
+    Projection    UTM
+    Zone          56
+    Datum         WGS84
+    Zunits        NO
+    Units         METERS
+    Spheroid      WGS84
+    Xshift        0.0000000000
+    Yshift        10000000.0000000000
+    Parameters
     """ 
 
@@ -945 +977 @@
 
     root, ext = os.path.splitext(filename)
-    fid = open(filename)
+
+    ###########################################
+    # Read Meta data
+    if verbose: print 'Reading METADATA from %s' %root + '.prj'
+    metadatafile = open(root + '.prj')    
+    metalines = metadatafile.readlines()
+    metadatafile.close()
+
+    L = metalines[0].strip().split()
+    assert L[0].strip().lower() == 'projection'
+    projection = L[1].strip()                   #TEXT
+
+    L = metalines[1].strip().split()
+    assert L[0].strip().lower() == 'zone'
+    zone = int(L[1].strip())    
+
+    L = metalines[2].strip().split()
+    assert L[0].strip().lower() == 'datum'
+    datum = L[1].strip()                        #TEXT 
+
+    L = metalines[3].strip().split()      
+    assert L[0].strip().lower() == 'zunits'     #IGNORE
+    zunits = L[1].strip()                       #TEXT
+
+    L = metalines[4].strip().split()
+    assert L[0].strip().lower() == 'units'
+    units = L[1].strip()                        #TEXT
+
+    L = metalines[5].strip().split()
+    assert L[0].strip().lower() == 'spheroid'   #IGNORE
+    spheroid = L[1].strip()                     #TEXT
+
+    L = metalines[6].strip().split()
+    assert L[0].strip().lower() == 'xshift'     
+    false_easting = float(L[1].strip())
+
+    L = metalines[7].strip().split()
+    assert L[0].strip().lower() == 'yshift'
+    false_northing = float(L[1].strip())    
+
+    #print false_easting, false_northing, zone, datum
+
+
+    ###########################################
+    #Read DEM data
+    datafile = open(filename)
 
     if verbose: print 'Reading DEM from %s' %filename
-    lines = fid.readlines()
-    fid.close()
+    lines = datafile.readlines()
+    datafile.close()
 
     if verbose: print 'Got', len(lines), ' lines'
@@ -962 +1039 @@
     assert len(lines) == nrows + 6
 
+
+    ##########################################
+
+    
     netcdfname = root + '.dem'
     if verbose: print 'Store to NetCDF file %s' %netcdfname
@@ -979 +1060 @@
     fid.NODATA_value = NODATA_value
 
+    fid.zone = zone
+    fid.false_easting = false_easting
+    fid.false_northing = false_northing    
+    fid.projection = projection
+    fid.datum = datum
+    fid.units = units
+    
+
     # dimension definitions
     fid.createDimension('number_of_rows', nrows)
@@ -1003 +1092 @@
 
 
-
-
-
-
 def ferret2sww(basefilename, verbose=False,
                minlat = None, maxlat =None,
                minlon = None, maxlon =None,
                mint = None, maxt = None, mean_stage = 0,
-               origin = (0,0,0)):
+               origin = (0,0,0,0,0)):
     """Convert 'Ferret' NetCDF format for wave propagation to
     sww format native to pyvolution.
@@ -1169 +1254 @@
 
     #Check zone boundaries
-    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
+    refzone, _, _ = redfearn(latitudes[0],longitudes[0],
+                             false_easting=origin[3],
+                             false_northing=origin[4])
+    
+                             
 
     vertices = {}
@@ -1176 +1265 @@
             vertices[l,k] = i
             
-            zone, easting, northing = redfearn(lat,lon)
+            zone, easting, northing = redfearn(lat,lon,
+                                               false_easting=origin[3],
+                                               false_northing=origin[4])
 
             msg = 'Zone boundary crossed at longitude =', lon