Package SatStress :: Module GridCalc
[frames] | no frames]

Source Code for Module SatStress.GridCalc

  1  """ 
  2  Calculate stresses on an icy satellite over a rectangular geographic region on 
  3  a regularly spaced lat-lon grid. 
  4   
  5  The datacube containing the results of the calculation are output as a 
  6  U{Unidata netCDF <http://www.unidata.ucar.edu/software/netcdf>} (.nc) file, 
  7  which can be displayed using a wide variety of visualization software. 
  8   
  9  """ 
 10   
 11  import SatStress as SS 
 12  import time 
 13  import netCDF3 
 14  import physcon as pc 
 15  import numpy 
 16  from numpy import linalg 
 17  import scipy 
 18  from optparse import OptionParser 
 19   
 20   
21 -def main():
22 """ 23 Calculate satellite stresses on a regular grid within a lat-lon window. 24 25 """ 26 27 usage = "usage: %prog [options] satfile gridfile outfile" 28 description = __doc__ 29 op = OptionParser(usage) 30 31 (options, args) = op.parse_args() 32 33 # Do a little error checking on the arguments: 34 if len(args) != 3: 35 op.error("incorrect number of arguments") 36 37 # The satfile defines the satellite and forcings it is subject to: 38 satfile = open(args[0], 'r') 39 the_sat = SS.Satellite(satfile) 40 satfile.close() 41 42 # The gridfile defines the scope and resolution of the calculation: 43 gridfile = open(args[1], 'r') 44 the_grid = Grid(gridfile, satellite=the_sat) 45 gridfile.close() 46 47 # Initialize a StressCalc object for the diurnal and NSR stresses 48 diurnal_stress = SS.StressCalc([SS.Diurnal(the_sat),]) 49 nsr_stress = SS.StressCalc([SS.NSR(the_sat),]) 50 51 # Make a netCDF object, and set its metadata 52 # Create a netCDF file object to stick the calculation results in: 53 nc_out = netCDF3.Dataset(args[2], 'w') 54 55 # Set metadata fields of nc_out appropriate to the calculation at hand. 56 57 # - Make sure that the file complies with either COARDS or CF. 58 # - Include all of the parameters that go into doing the calculation (whatever 59 # is listed within the .satellite and .grid files). 60 # - Make sure that the units for those parameters are clear. 61 # - A decent description of the contents of the file in English. 62 63 # Set the global data attributes of the netCDF file describing the run 64 # This allows us to re-construct the entire model run if we want to, given 65 # the SatStress library and a netCDF output file. 66 nc_out.description = "Testing pySatStress on a regular grid" 67 nc_out.history = "Created: %s using pySatStress" % ( time.ctime(time.time()) ) 68 nc_out.Conventions = "COARDS" 69 70 # Specify the size and shape of the output datacube 71 72 n_deltas = the_grid.nsr_delta_numsteps 73 nc_out.createDimension('nsr_delta', n_deltas) 74 nsr_deltas = nc_out.createVariable('nsr_delta', 'f4', ('nsr_delta',)) 75 nsr_deltas.units = "" 76 nsr_deltas.long_name = "NSR Surface Delta (mu/(eta*omega))" 77 nsr_deltas[:] = numpy.logspace(the_grid.nsr_delta_min, the_grid.nsr_delta_max, num=n_deltas) 78 79 n_times = len(numpy.arange(the_grid.time_min, the_grid.time_max, the_grid.time_step)) 80 nc_out.createDimension('time', n_times) 81 times = nc_out.createVariable('time', 'f4', ('time',)) 82 83 # Check to see what kind of units we're using for time, 84 # and name the variables and their units appropriately 85 if the_grid.orbit_min is None: 86 times.units = "seconds" 87 times.long_name = "time after periapse" 88 times[:] = numpy.arange(the_grid.time_min, the_grid.time_max, the_grid.time_step) 89 else: 90 times.units = "degrees" 91 times.long_name = "degrees after periapse" 92 times[:] = numpy.arange(the_grid.orbit_min, the_grid.orbit_max, the_grid.orbit_step) 93 94 n_lats = len(numpy.arange(the_grid.lat_min, the_grid.lat_max+the_grid.latlon_step, the_grid.latlon_step)) 95 nc_out.createDimension('latitude', n_lats) 96 lats = nc_out.createVariable('latitude', 'f4', ('latitude',)) 97 lats.units = "degrees_north" 98 lats.long_name = "latitude" 99 lats[:] = numpy.arange(the_grid.lat_min, the_grid.lat_max + (the_grid.latlon_step/2.0), the_grid.latlon_step) 100 101 n_lons = len(numpy.arange(the_grid.lon_min, the_grid.lon_max+the_grid.latlon_step, the_grid.latlon_step)) 102 nc_out.createDimension('longitude', n_lons) 103 lons = nc_out.createVariable('longitude', 'f4', ('longitude',)) 104 lons.units = "degrees_east" 105 lons.long_name = "longitude" 106 lons[:] = numpy.arange(the_grid.lon_min, the_grid.lon_max + (the_grid.latlon_step/2.0), the_grid.latlon_step) 107 # end dimension and coordinate variable definitions 108 109 # Define the data variables 110 Ttt_D = nc_out.createVariable('Ttt_D', 'f4', ('time', 'latitude', 'longitude',)) 111 Ttt_D.units = "Pa" 112 Ttt_D.long_name = "north-south component stress of Diurnal stresses" 113 114 Tpt_D = nc_out.createVariable('Tpt_D', 'f4', ('time', 'latitude', 'longitude',)) 115 Tpt_D.units = "Pa" 116 Tpt_D.long_name = "shear component of Diurnal stresses" 117 118 Tpp_D = nc_out.createVariable('Tpp_D', 'f4', ('time', 'latitude', 'longitude',)) 119 Tpp_D.units = "Pa" 120 Tpp_D.long_name = "east-west component of Diurnal stresses" 121 122 Ttt_N = nc_out.createVariable('Ttt_N', 'f4', ('nsr_delta', 'latitude', 'longitude',)) 123 Ttt_N.units = "Pa" 124 Ttt_N.long_name = "north-south component of NSR stresses" 125 126 Tpt_N = nc_out.createVariable('Tpt_N', 'f4', ('nsr_delta', 'latitude', 'longitude',)) 127 Tpt_N.units = "Pa" 128 Tpt_N.long_name = "shear component of NSR stresses" 129 130 Tpp_N = nc_out.createVariable('Tpp_N', 'f4', ('nsr_delta', 'latitude', 'longitude',)) 131 Tpp_N.units = "Pa" 132 Tpp_N.long_name = "east-west component of NSR stresses" 133 # end data variable definitions 134 135 # Loop over the time variable, doing diurnal calculations over an orbit 136 for t in range(n_times): 137 # We need some kind of progress update, and we need to make sure that 138 # we have a representation of the time coordinate in seconds, because 139 # that's what the SatStress library expects - even if we're ultimately 140 # communicating time to the user in terms of "degrees after periapse" 141 if the_grid.orbit_min is None: 142 time_sec = times[t] 143 else: 144 time_sec = the_sat.orbit_period()*(times[t]/360.0) 145 146 print "Calculating diurnal stresses at", times[t], times.long_name 147 for lon in range(n_lons): 148 for lat in range(n_lats): 149 150 Tau_D = diurnal_stress.tensor(theta = scipy.radians(90.0-lats[lat]),\ 151 phi = scipy.radians(lons[lon]),\ 152 t = time_sec ) 153 154 nc_out.variables['Ttt_D'][t,lat,lon] = Tau_D[0,0] 155 nc_out.variables['Tpt_D'][t,lat,lon] = Tau_D[1,0] 156 nc_out.variables['Tpp_D'][t,lat,lon] = Tau_D[1,1] 157 # Not doing the principal components in this app... but save for later. 158 #(magA, magB), eigenvecs = linalg.eigh(Tau_D) 159 #mag1_D[t,lat,lon], mag2_D[t,lat,lon] = max(magA, magB), min(magA, magB) 160 161 # Make sure everything gets written out to the file. 162 nc_out.sync() 163 # end Diurnal calculations 164 165 # Begin the NSR calculations: 166 for delta in range(n_deltas): 167 # Set the nsr_period such that it corresponds to the current value of Delta: 168 the_sat.nsr_period = 4*numpy.pi*nsr_deltas[delta]*the_sat.layers[-1].viscosity/the_sat.layers[-1].lame_mu 169 # Re-set the StressCalc object based on the new nsr_period: 170 nsr_stress = SS.StressCalc([SS.NSR(the_sat),]) 171 # We need some kind of progress update... 172 print "Calculating NSR stresses at Delta =", nsr_deltas[delta], "(P_nsr = %g yr)" % (the_sat.nsr_period/31556926) 173 for lon in range(n_lons): 174 for lat in range(n_lats): 175 176 # Note that t is always zero, because in the NSR potential, all time 177 # does is translate surface features through the field, and here we're 178 # actually not interested in what the surface features are doing - we 179 # just want to know what the stresses look like relative to a fixed 180 # longitude (relative to the planet's location in the satellite's sky) 181 # There's more discussion of this issue in Wahr et al. (2008). 182 Tau_N = nsr_stress.tensor(theta = scipy.radians(90.0-lats[lat]),\ 183 phi = scipy.radians(lons[lon]),\ 184 t = 0) 185 186 nc_out.variables['Ttt_N'][delta,lat,lon] = Tau_N[0,0] 187 nc_out.variables['Tpt_N'][delta,lat,lon] = Tau_N[1,0] 188 nc_out.variables['Tpp_N'][delta,lat,lon] = Tau_N[1,1] 189 190 # Not doing the principal components in this app... but save for later. 191 #(magA, magB), eigenvecs = linalg.eigh(Tau_N) 192 #mag1_N[delta,lat,lon], mag2_N[delta,lat,lon] = max(magA, magB), min(magA, magB) 193 194 # re-set the NSR period to correspond to the next value of Delta... 195 196 # end NSR calculation loop. 197 198 # Make sure everything gets written out to the file. 199 nc_out.sync()
200 201 # end main() 202
203 -class Grid(): #
204 """ 205 A container class defining the temporal and geographic range and resolution 206 of the calculation. 207 208 The parameters defining the calculation grid are read in from a name value 209 file, parsed into a Python dictionary using L{SatStress.nvf2dict}, and used 210 to set the data attributes of the L{Grid} object. 211 212 The geographic extent of the calculation is specified by minimum and maximum 213 values for latitude and longitude. 214 215 The geographic resolution of the calculation is defined by an angular 216 separation between calculations. This angular separation is the same in 217 the north-south and the east-west direction. 218 219 The temporal range and resolution of the calculation can be specified either 220 in terms of actual time units (seconds) or in terms of the satellite's orbital 221 position (in degrees). In both cases, time=0 is taken to occur at periapse. 222 223 S{Delta} is a measure of how viscous or elastic the response of the body 224 is. It's equal to (S{mu})/(S{eta}S{omega}) where S{mu} and S{eta} are the 225 shear modulus and viscosity of the surface layer, respectively, and 226 S{omega} is the forcing frequency to which the body is subjected (see Wahr 227 et al. (2008) for a detailed discussion). It is a logarithmic parameter, so 228 its bounds are specified as powers of 10, e.g. if the minimum value is -3, 229 the initial S{Delta} is 10^-3 = 0.001. 230 231 @ivar grid_id: A string identifying the grid 232 @type grid_id: str 233 @ivar lat_min: Southern bound, degrees (north positive). 234 @type lat_min: float 235 @ivar lat_max: Northern bound, degrees (north positive). 236 @type lat_max: float 237 @ivar lon_min: Western bound, degrees (east positive). 238 @type lon_min: float 239 @ivar lon_max: Eastern bound, degrees (east positive). 240 @type lon_max: float 241 @ivar latlon_step: Angular separation between calculations. 242 @type latlon_step: float 243 @ivar time_min: Initial time at which calculation begins (0 = periapse). 244 @type time_min: float 245 @ivar time_max: Final time at which calculation ends. 246 @type time_max: float 247 @ivar time_step: Seconds between subsequent calculations. 248 @type time_step: float 249 @ivar orbit_min: Initial orbital position in degrees (0 = periapse) 250 @type orbit_min: float 251 @ivar orbit_max: Final orbital position in degrees (0 = periapse) 252 @type orbit_max: float 253 @ivar orbit_step: Orbital angular separation between calculations in degrees 254 @type orbit_step: float 255 @ivar nsr_delta_min: Initial S{Delta} = 10^(nsr_delta_min) 256 @type nsr_delta_min: float 257 @ivar nsr_delta_max: Final S{Delta} = 10^(nsr_delta_max) 258 @type nsr_delta_max: float 259 @ivar nsr_delta_numsteps: How many S{Delta} values to calculate total 260 @type nsr_delta_numsteps: int 261 262 """ 263
264 - def __init__(self, gridFile, satellite=None):
265 """Initialize the Grid object from a gridFile.""" 266 267 # Slurp up the proffered parameters: 268 gridParams = SS.nvf2dict(gridFile, comment='#') 269 270 self.grid_id = gridParams['GRID_ID'] 271 272 self.lat_min = float(gridParams['LAT_MIN']) 273 self.lat_max = float(gridParams['LAT_MAX']) 274 self.lon_min = float(gridParams['LON_MIN']) 275 self.lon_max = float(gridParams['LON_MAX']) 276 self.latlon_step = float(gridParams['LATLON_STEP']) 277 278 self.nsr_delta_numsteps = float(gridParams['NSR_DELTA_NUMSTEPS']) 279 self.nsr_delta_min = float(gridParams['NSR_DELTA_MIN']) 280 self.nsr_delta_max = float(gridParams['NSR_DELTA_MAX']) 281 282 # We MUST have all three of either of these and NONE of the other. 283 if gridParams.has_key('TIME_MIN'): 284 self.time_min = float(gridParams['TIME_MIN']) 285 self.time_max = float(gridParams['TIME_MAX']) 286 self.time_step = float(gridParams['TIME_STEP']) 287 288 self.orbit_min = None 289 self.orbit_max = None 290 self.orbit_step = None 291 292 elif gridParams.has_key('ORBIT_MIN'): 293 self.orbit_min = float(gridParams['ORBIT_MIN']) 294 self.orbit_max = float(gridParams['ORBIT_MAX']) 295 self.orbit_step = float(gridParams['ORBIT_STEP']) 296 297 self.time_min = satellite.orbit_period()*(self.orbit_min/360.0) 298 self.time_max = satellite.orbit_period()*(self.orbit_max/360.0) 299 self.time_step = satellite.orbit_period()*(self.orbit_step/360.0) 300 301 else: 302 # Die with an error about bad input file... 303 pass
304 305 # end class Grid 306 307 if __name__ == "__main__": 308 main() 309