tPISMNC.py - pism-exp-gsw - ice stream and sediment transport experiments
 (HTM) git clone git://src.adamsgaard.dk/pism-exp-gsw
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       tPISMNC.py (5381B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 try:
            4     import netCDF4 as netCDF
            5 except:
            6     print("netCDF4 is not installed!")
            7     sys.exit(1)
            8 
            9 
           10 class PISMDataset(netCDF.Dataset):
           11 
           12     def create_time(self, use_bounds=False, length=None, units=None):
           13         self.createDimension('time', size=length)
           14         t_var = self.createVariable('time', 'f8', ('time',))
           15 
           16         t_var.axis = "T"
           17         t_var.long_name = "time"
           18         if not units:
           19             t_var.units = "seconds since 1-1-1"  # just a default
           20         else:
           21             t_var.units = units
           22 
           23         if use_bounds:
           24             self.createDimension('n_bounds', 2)
           25             self.createVariable("time_bounds", 'f8', ('time', 'n_bounds'))
           26             t_var.bounds = "time_bounds"
           27 
           28     def create_dimensions(self, x, y, time_dependent=False, use_time_bounds=False):
           29         """
           30         Create PISM-compatible dimensions in a NetCDF file.
           31         """
           32 
           33         if time_dependent and not 'time' in list(self.variables.keys()):
           34             self.create_time(use_time_bounds)
           35 
           36         self.createDimension('x', x.size)
           37         self.createDimension('y', y.size)
           38 
           39         x_var = self.createVariable('x', 'f8', ('x',))
           40         x_var[:] = x
           41 
           42         y_var = self.createVariable('y', 'f8', ('y',))
           43         y_var[:] = y
           44 
           45         x_var.axis = "X"
           46         x_var.long_name = "X-coordinate in Cartesian system"
           47         x_var.units = "m"
           48         x_var.standard_name = "projection_x_coordinate"
           49 
           50         y_var.axis = "Y"
           51         y_var.long_name = "Y-coordinate in Cartesian system"
           52         y_var.units = "m"
           53         y_var.standard_name = "projection_y_coordinate"
           54 
           55         self.sync()
           56 
           57     def append_time(self, value, bounds=None):
           58         if 'time' in list(self.dimensions.keys()):
           59             time = self.variables['time']
           60             N = time.size
           61             time[N] = value
           62 
           63             if bounds:
           64                 self.variables['time_bounds'][N, :] = bounds
           65 
           66     def set_attrs(self, var_name, attrs):
           67         """attrs should be a list of (name, value) tuples."""
           68         if not attrs:
           69             return
           70 
           71         for (name, value) in attrs.items():
           72             if name == "_FillValue":
           73                 continue
           74             setattr(self.variables[var_name], name, value)
           75 
           76     def define_2d_field(self, var_name, time_dependent=False, dims=None, nc_type='f8', attrs=None):
           77         """
           78         time_dependent: boolean
           79 
           80         dims: an optional list of dimension names. use this to override the
           81               default order ('time', 'y', 'x')
           82 
           83         attrs: a dictionary of attributes
           84         """
           85         if not dims:
           86             if time_dependent:
           87                 dims = ('time', 'y', 'x')
           88             else:
           89                 dims = ('y', 'x')
           90 
           91         try:
           92             var = self.variables[var_name]
           93         except:
           94             if attrs is not None and '_FillValue' in list(attrs.keys()):
           95                 var = self.createVariable(var_name, nc_type, dims,
           96                                           fill_value=attrs['_FillValue'])
           97             else:
           98                 var = self.createVariable(var_name, nc_type, dims)
           99 
          100             self.set_attrs(var_name, attrs)
          101 
          102         return var
          103 
          104     def define_timeseries(self, var_name, attrs=None):
          105         try:
          106             if attrs is not None and '_FillValue' in list(attrs.keys()):
          107                 var = self.createVariable(var_name, 'f8', ('time',),
          108                                           fill_value=attrs['_FillValue'])
          109             else:
          110                 var = self.createVariable(var_name, 'f8', ('time',))
          111         except:
          112             var = self.variables[var_name]
          113 
          114         self.set_attrs(var_name, attrs)
          115 
          116         return var
          117 
          118     def write(self, var_name, data, time_dependent=False, attrs=None):
          119         """
          120         Write time-series or a 2D field to a file.
          121         """
          122 
          123         if data.ndim == 1:
          124             return self.write_timeseries(var_name, data, attrs=attrs)
          125         elif data.ndim == 2:
          126             return self.write_2d_field(var_name, data, time_dependent, attrs=attrs)
          127         else:
          128             return None
          129 
          130     def write_2d_field(self, var_name, data, time_dependent=False, attrs=None):
          131         """
          132         Write a 2D numpy array to a file in a format PISM can read.
          133         """
          134 
          135         var = self.define_2d_field(var_name, time_dependent, attrs=attrs)
          136 
          137         if time_dependent:
          138             last_record = self.variables['time'].size - 1
          139             var[last_record, :, :] = data
          140         else:
          141             var[:] = data
          142 
          143         return var
          144 
          145     def write_timeseries(self, var_name, data, attrs=None):
          146         """Write a 1D (time-series) array to a file."""
          147 
          148         var = self.define_timeseries(var_name, attrs=attrs)
          149         var[:] = data
          150 
          151         return var
          152 
          153 
          154 if __name__ == "__main__":
          155     # produce a NetCDF file for testing
          156     from numpy import linspace, meshgrid
          157 
          158     nc = PISMDataset("foo.nc", 'w')
          159 
          160     x = linspace(-100, 100, 101)
          161     y = linspace(-100, 100, 201)
          162 
          163     xx, yy = meshgrid(x, y)
          164 
          165     nc.create_dimensions(x, y, time_dependent=True, use_time_bounds=True)
          166 
          167     nc.define_2d_field("xx", time_dependent=True,
          168                        attrs={"long_name": "xx",
          169                               "comment": "test variable",
          170                               "valid_range": (-200.0, 200.0)})
          171 
          172     for t in [0, 1, 2, 3]:
          173         nc.append_time(t, (t - 1, t))
          174 
          175         nc.write("xx", xx + t, time_dependent=True)
          176         nc.write("yy", yy + 2 * t, time_dependent=True)
          177 
          178     nc.close()