Changeset 4387

Timestamp:

Apr 17, 2007, 5:07:33 PM (18 years ago)

Author:

duncan

Message:

small modifications to geo-reffing

Location:

anuga_core/source/anuga

Files:

• coordinate_transforms/geo_reference.py (modified) (7 diffs)
• shallow_water/data_manager.py (modified) (10 diffs)
• shallow_water/test_data_manager.py (modified) (2 diffs)

anuga_core/source/anuga/coordinate_transforms/geo_reference.py

@@ -16 +16 @@
 TITLE = '#geo reference' + "\n" #this title is referred to in the .xya format
 
+DEFAULT_PROJECTION = 'UTM'
+DEFAULT_DATUM = 'wgs84'
+DEFAULT_UNITS = 'm'
+DEFAULT_FALSE_EASTING = 500000
+DEFAULT_FALSE_NORTHING = 10000000 #Default for southern hemisphere 
+
 class Geo_reference:
     """
@@ -24 +30 @@
                  xllcorner = 0.0,
                  yllcorner = 0.0,
-                 datum = 'wgs84',
-                 projection = 'UTM',                 units = 'm',
-                 false_easting = 500000,
-                 false_northing = 10000000, #Default for southern hemisphere 
+                 datum = DEFAULT_DATUM,
+                 projection = DEFAULT_PROJECTION,
+                 units = DEFAULT_UNITS,
+                 false_easting = DEFAULT_FALSE_EASTING,
+                 false_northing = DEFAULT_FALSE_NORTHING, 
                  NetCDFObject=None,
                  ASCIIFile=None,
@@ -39 +46 @@
          it can't unread it, so this info has to be passed.
          If you know of a way to unread this info, then tell us.
-        """
+
+         Note, the text file only saves a sub set of the info the
+         points file does.  Currently the info not written in text
+         must be the default info, since ANUGA assumes it isn't
+         changing.
+         """
 
         self.false_easting = false_easting
@@ -82 +94 @@
         self.zone = infile.zone[0]
 
+        self.false_easting = infile.false_easting[0]
+        self.false_northing = infile.false_northing[0]
+        
+        self.datum = infile.datum        
+        self.projection = infile.projection
+        self.units = infile.units
+        
         # Fix some assertion failures
         if type(self.zone) == ArrayType and self.zone.shape == ():
@@ -95 +114 @@
                 type(self.yllcorner) == types.IntType)
         assert (type(self.zone) == types.IntType)
-
-        #FIXME (Ole) Read in the rest...
+        
+        assert (self.false_easting == DEFAULT_FALSE_EASTING)
+        assert (self.false_northing == DEFAULT_FALSE_NORTHING)
+
+        assert(self.datum == DEFAULT_DATUM)
+        assert(self.projection == DEFAULT_PROJECTION)
+        assert (self.units == DEFAULT_UNITS)
         
         
@@ -132 +156 @@
         assert (type(self.zone) == types.IntType)
         
-        #false_easting = infile.false_easting[0]
-        #false_northing = infile.false_northing[0]
         
     def change_points_geo_ref(self, points,
@@ -267 +289 @@
             cmp = 1
         return cmp
-        
+
+def write_NetCDF_georeference(origin, outfile):
+    """
+    Write georeferrence info to a netcdf file, usually sww.
+
+    The origin can be a georef instance or parrameters for a geo_ref instance
+
+    outfile is the name of the file to be written to.
+    """
+    if isinstance(origin, Geo_reference): 
+        geo_ref = origin
+    else:
+        geo_ref = apply(Geo_reference,origin)
+    geo_ref.write_NetCDF(outfile)
+    return geo_ref
+    
 #-----------------------------------------------------------------------
 

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -72 +72 @@
 from anuga.coordinate_transforms.redfearn import redfearn, \
      convert_from_latlon_to_utm
-from anuga.coordinate_transforms.geo_reference import Geo_reference
+from anuga.coordinate_transforms.geo_reference import Geo_reference, \
+     write_NetCDF_georeference
 from anuga.geospatial_data.geospatial_data import Geospatial_data,\
      ensure_absolute
@@ -2645 +2646 @@
     outfile = NetCDFFile(swwname, 'w')
 
-    #Create new file
-    outfile.institution = 'Geoscience Australia'
-    outfile.description = 'Converted from Ferret files: %s, %s, %s, %s'\
-                          %(basename_in + '_ha.nc',
-                            basename_in + '_ua.nc',
-                            basename_in + '_va.nc',
-                            basename_in + '_e.nc')
-
-
-    #For sww compatibility
-    outfile.smoothing = 'Yes'
-    outfile.order = 1
-
-    #Start time in seconds since the epoch (midnight 1/1/1970)
-    outfile.starttime = starttime = times[0]
-    times = times - starttime  #Store relative times
-
-    # dimension definitions
-    outfile.createDimension('number_of_volumes', number_of_volumes)
-
-    outfile.createDimension('number_of_vertices', 3)
-    outfile.createDimension('number_of_points', number_of_points)
-
-
-    #outfile.createDimension('number_of_timesteps', len(times))
-    outfile.createDimension('number_of_timesteps', len(times))
-
-    # variable definitions
-    outfile.createVariable('x', precision, ('number_of_points',))
-    outfile.createVariable('y', precision, ('number_of_points',))
-    outfile.createVariable('elevation', precision, ('number_of_points',))
-
-    #FIXME: Backwards compatibility
-    outfile.createVariable('z', precision, ('number_of_points',))
-    #################################
-
-    outfile.createVariable('volumes', Int, ('number_of_volumes',
-                                            'number_of_vertices'))
-
-    outfile.createVariable('time', precision,
-                           ('number_of_timesteps',))
-
-    outfile.createVariable('stage', precision,
-                           ('number_of_timesteps',
-                            'number_of_points'))
-
-    outfile.createVariable('xmomentum', precision,
-                           ('number_of_timesteps',
-                            'number_of_points'))
-
-    outfile.createVariable('ymomentum', precision,
-                           ('number_of_timesteps',
-                            'number_of_points'))
-
+    description = 'Converted from Ferret files: %s, %s, %s, %s'\
+                  %(basename_in + '_ha.nc',
+                    basename_in + '_ua.nc',
+                    basename_in + '_va.nc',
+                    basename_in + '_e.nc')
+    write_sww_header(outfile, times, number_of_volumes,
+                     number_of_points, description=description) 
+#     #Create new file
+#     outfile.institution = 'Geoscience Australia'
+#     outfile.description = 'Converted from Ferret files: %s, %s, %s, %s'\
+#                           %(basename_in + '_ha.nc',
+#                             basename_in + '_ua.nc',
+#                             basename_in + '_va.nc',
+#                             basename_in + '_e.nc')
+
+
+#     #For sww compatibility
+#     outfile.smoothing = 'Yes'
+#     outfile.order = 1
+
+#     #Start time in seconds since the epoch (midnight 1/1/1970)
+#     outfile.starttime = starttime = times[0]
+#     times = times - starttime  #Store relative times
+
+#     # dimension definitions
+#     outfile.createDimension('number_of_volumes', number_of_volumes)
+
+#     outfile.createDimension('number_of_vertices', 3)
+#     outfile.createDimension('number_of_points', number_of_points)
+
+
+#     #outfile.createDimension('number_of_timesteps', len(times))
+#     outfile.createDimension('number_of_timesteps', len(times))
+
+#     # variable definitions
+#     outfile.createVariable('x', precision, ('number_of_points',))
+#     outfile.createVariable('y', precision, ('number_of_points',))
+#     outfile.createVariable('elevation', precision, ('number_of_points',))
+
+#     #FIXME: Backwards compatibility
+#     outfile.createVariable('z', precision, ('number_of_points',))
+#     #################################
+
+#     outfile.createVariable('volumes', Int, ('number_of_volumes',
+#                                             'number_of_vertices'))
+
+#     outfile.createVariable('time', precision,
+#                            ('number_of_timesteps',))
+
+#     outfile.createVariable('stage', precision,
+#                            ('number_of_timesteps',
+#                             'number_of_points'))
+
+#     outfile.createVariable('xmomentum', precision,
+#                            ('number_of_timesteps',
+#                             'number_of_points'))
+
+#     outfile.createVariable('ymomentum', precision,
+#                            ('number_of_timesteps',
+#                             'number_of_points'))
+
+#     outfile.variables['time'][:] = times   #Store time relative
 
     #Store
@@ -2741 +2750 @@
     volumes = array(volumes)
 
-    if origin == None:
-        zone = refzone
-        xllcorner = min(x)
-        yllcorner = min(y)
-    else:
-        zone = origin[0]
-        xllcorner = origin[1]
-        yllcorner = origin[2]
-
-
-    outfile.xllcorner = xllcorner
-    outfile.yllcorner = yllcorner
-    outfile.zone = zone
+    if origin is None:
+        origin = Geo_reference(refzone,min(x),min(y))
+    geo_ref = write_NetCDF_georeference(origin, outfile)
+    
+#     if origin == None:
+#         zone = refzone
+#         xllcorner = min(x)
+#         yllcorner = min(y)
+#     else:
+#         zone = origin[0]
+#         xllcorner = origin[1]
+#         yllcorner = origin[2]
+
+
+#     outfile.xllcorner = xllcorner
+#     outfile.yllcorner = yllcorner
+#     outfile.zone = zone
 
 
@@ -2767 +2780 @@
     from Numeric import resize
     z = resize(z,outfile.variables['z'][:].shape)
-    outfile.variables['x'][:] = x - xllcorner
-    outfile.variables['y'][:] = y - yllcorner
+    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
+    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
     outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
     outfile.variables['elevation'][:] = z
-    outfile.variables['time'][:] = times   #Store time relative
     outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
 
@@ -2808 +2820 @@
         print '  Reference:'
         print '    Lower left corner: [%f, %f]'\
-              %(xllcorner, yllcorner)
+              %(geo_ref.get_xllcorner(), geo_ref.get_yllcorner())
         print '    Start time: %f' %starttime
         print '  Extent:'
@@ -4573 +4585 @@
     return mux_times_start_i, mux_times_fin_i
 
-def write_sww_header(outfile, times, number_of_volumes, number_of_points ):
+def write_sww_header(outfile, times, number_of_volumes, number_of_points, description='Converted from XXX'):
     """
     outfile - the name of the file that will be written
@@ -4579 +4591 @@
     number_of_volumes - the number of triangles
     """
+    times = ensure_numeric(times)
+    
     number_of_times = len(times)
     
     outfile.institution = 'Geoscience Australia'
     outfile.description = 'Converted from XXX'
-
 
     #For sww compatibility
@@ -4594 +4607 @@
     else:
         outfile.starttime = starttime = times[0]
-
+    times = times - starttime  #Store relative times
 
     # dimension definitions
     outfile.createDimension('number_of_volumes', number_of_volumes)
-
     outfile.createDimension('number_of_vertices', 3)
     outfile.createDimension('number_of_points', number_of_points)
@@ -4630 +4642 @@
                             'number_of_points'))
     
-    outfile.variables['time'][:] = times   
+    outfile.variables['time'][:] = times    #Store time relative
 
 
@@ -4639 +4651 @@
     volumes = array(volumes)
 
-    if origin is not None:
-        if isinstance(origin, Geo_reference): 
-            geo_ref = origin
-        else:
-            geo_ref = apply(Geo_reference,origin)
-    else:
-        geo_ref = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1]))
-    #geo_ref = Geo_reference(zone,0,0)
-    geo_ref.write_NetCDF(outfile)
+    if origin is None:
+        origin = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1]))
+    geo_ref = write_NetCDF_georeference(origin, outfile)
+    
+#     if origin is not None:
+#         if isinstance(origin, Geo_reference): 
+#             geo_ref = origin
+#         else:
+#             geo_ref = apply(Geo_reference,origin)
+#     else:
+#         geo_ref = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1]))
+#     #geo_ref = Geo_reference(zone,0,0)
+#     geo_ref.write_NetCDF(outfile)
+    
 
     # This will put the geo ref in the middle

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -5724 +5724 @@
         times = fid.variables['time'][:]
         
-        times_actual = [200,300,400,500]
+        times_actual = [0,100,200,300]
        
         assert allclose(ensure_numeric(times),
@@ -5974 +5974 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww_origin')
     #suite = unittest.makeSuite(Test_Data_Manager,'cache_test_URS_points_needed_and_urs_ungridded2sww')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs_ungridded_hole')