5.4. Make 4D fields (to netCDF)

5.4.1. with integrated functionalities: meta-resources

In [30]:
import os
import epygram
epygram.init_env()
workdir = epygram.config.userlocaldir + '/notebooks_data'
os.chdir(workdir)
import numpy
In [31]:
files = ['advanced_examples/ICMSHAROM+{:0>4}'.format(i) for i in range(0,2+1)]
In [32]:
# meta-resource are a wrapper over raw resources, enabling them to combine H2D fields on a vertical discretization
# as 3D fields, or combine fields from a temporal series of resources as a field with a temporal dimension
r = epygram.resources.meta_resource(files, 'r', 'MV+CL')  # MV = Multi-Validities, CL = Combine-Levels

# such meta-resources (CL) need a 'generic' or GRIB2 kind of fid to read in, for recombining coherent levels
# of course it underneath needs a fid conversion, which is implemented basically only for FA and LFI
t = r.readfield({'discipline':0, 'parameterCategory':0, 'parameterNumber':0,  # temperature, as GRIB2
                 'typeOfFirstFixedSurface':119})  # on hybrid-pressure (model) levels
In [33]:
print(t.data.shape)
t.sp2gp()
print(t.data.shape)
(3, 60, 423256)
(3, 60, 720, 750)
In [34]:
rh = r.readfield({'discipline':0, 'parameterCategory':1, 'parameterNumber':1,
                  'typeOfFirstFixedSurface':103, 'level':2})  # relative humidity at 2m, as GRIB2
In [35]:
out = epygram.formats.resource('out4D.nc', 'w', fmt='netCDF')  # open the output netCDF
In [36]:
t.fid['netCDF'] = 'temperature'
out.writefield(t)
In [37]:
rh.fid['netCDF'] = 'relative_humidity'
out.writefield(rh)
# [2017/08/25-10:52:38][epygram.formats.netCDF][writefield:1039][INFO]: assume lons/lats match.
# [2017/08/25-10:52:38][epygram.formats.netCDF][writefield:1144][INFO]: assume projection parameters match.
In [38]:
out.close()
out.open(openmode='r')
out.what()
### FORMAT: netCDF

NetCDF Global Attributes:
        made_with: u'epygram-1.2.10'
        Conventions: u'CF-1.6'
NetCDF dimension information:
        Name: time
                size: 3
        Name: Z
                size: 60
        Name: Y
                size: 720
        Name: X
                size: 750
        Name: Z+1
                size: 61
NetCDF variable information:
        Name: time
                dimensions: (u'time',)
                size: 3
                type: dtype('float64')
                units: u'seconds since 2015-02-16 00:00:00'
        Name: Z
                dimensions: ()
                size: 1.0
                type: dtype('int64')
                standard_name: u'atmosphere_hybrid_sigma_pressure_coordinate'
                positive: u'down'
                formula_terms: u'ap: hybrid_coef_A b: hybrid_coef_B ps: surface_air_pressure'
                short_name: u'hybrid-pressure'
        Name: hybrid_coef_A
                dimensions: (u'Z+1',)
                size: 61
                type: dtype('float64')
        Name: hybrid_coef_B
                dimensions: (u'Z+1',)
                size: 61
                type: dtype('float64')
        Name: longitude
                dimensions: (u'Y', u'X')
                size: 540000
                type: dtype('float64')
                units: u'degrees'
        Name: latitude
                dimensions: (u'Y', u'X')
                size: 540000
                type: dtype('float64')
                units: u'degrees'
        Name: Projection_parameters
                dimensions: ()
                size: 1.0
                type: dtype('int64')
                grid_mapping_name: u'lambert_conformal_conic'
                earth_radius: 6371229.0
                x_resolution: 2500.0
                y_resolution: 2500.0
                longitude_of_central_meridian: 2.0
                latitude_of_projection_origin: 45.800000000000004
                standard_parallel: 45.800000000000004
                false_easting: 922500.0
                false_northing: 885000.0000025304
        Name: temperature
                dimensions: (u'time', u'Z', u'Y', u'X')
                size: 97200000
                type: dtype('float64')
                grid_mapping: u'Projection_parameters'
                vertical_grid: u'Z'
        Name: relative_humidity
                dimensions: (u'time', u'Y', u'X')
                size: 1620000
                type: dtype('float64')
                grid_mapping: u'Projection_parameters'

5.4.2. or do it manually

In [39]:
files = ['advanced_examples/GRIDFRANGP0025r0_{:0>4}'.format(i) for i in range(0,2+1)]  # these are GRIB1 files
print(files)
['advanced_examples/GRIDFRANGP0025r0_0000', 'advanced_examples/GRIDFRANGP0025r0_0001', 'advanced_examples/GRIDFRANGP0025r0_0002']
In [40]:
# try with Combine-Levels meta-resource
r = epygram.resources.meta_resource(files[0], 'r', 'CL')
In [41]:
t = r.readfield({'indicatorOfParameter':11,
                 'indicatorOfTypeOfLevel':100})  # temperature on Pressure levels, as GRIB1

AssertionErrorTraceback (most recent call last)
<ipython-input-41-375c9ec9f5f3> in <module>()
      1 t = r.readfield({'indicatorOfParameter':11,
----> 2                  'indicatorOfTypeOfLevel':100})  # temperature on Pressure levels, as GRIB1

/home/mary/EPyGrAM/src/epygram/resources/CombineLevelsResource.pyc in readfield(self, handgrip, getdata)
    158         :param getdata: if False, do not read data but only metadata
    159         """
--> 160         result = self.readfields(handgrip=handgrip, getdata=getdata)
    161         if len(result) != 1:
    162             raise epygramError(str(len(result)) + "field(s) have been found, only one expected.")

/home/mary/EPyGrAM/src/epygram/resources/CombineLevelsResource.pyc in readfields(self, handgrip, getdata)
    171         """
    172         fieldset = FieldSet()
--> 173         cont = self._create_list()
    174         for fid in self.listfields(select=handgrip):
    175             found = None

/home/mary/EPyGrAM/src/epygram/resources/CombineLevelsResource.pyc in _create_list(self)
     97         for fid in self.resource.listfields(complete=True):
     98             assert 'generic' in fid, \
---> 99                    "Not able to combine levels if fields do not have 'generic' fids"
    100             original_fid = fid[fmtfid(self.resource.format, fid)]
    101             generic_fid = fid['generic']

AssertionError: Not able to combine levels if fields do not have 'generic' fids
In [42]:
# This did not worked because 'CL' meta-resource need generic/GRIB2 => let's do it manually
In [43]:
r = epygram.resources.meta_resource(files, 'r', 'MV')  # this works anyway
In [44]:
t = r.readfield({'indicatorOfParameter':11,
                 'indicatorOfTypeOfLevel':100,
                 'level':850})  # temperature, as GRIB1
In [45]:
t.data.shape
Out[45]:
(3, 601, 801)
In [46]:
# either, we could have aggregated time dimension manually:
temp_time = [epygram.formats.resource(f, 'r').readfield({'indicatorOfParameter':11,
                                                         'indicatorOfTypeOfLevel':100,
                                                         'level':850})
             for f in files]
fld_t = temp_time[0]
for fld in temp_time[1:]:
    fld_t.extend(fld)
print(fld_t.data.shape)
(3, 601, 801)
In [47]:
# now let's gather vertical levels
levels = r.listfields(select={'indicatorOfParameter':11,
                              'indicatorOfTypeOfLevel':100},
                      onlykey='level')
print(levels)
[1000, 500, 600, 700, 800, 100, 150, 850, 200, 250, 900, 925, 950, 300, 400]
In [48]:
temp_levels = []  # we can read levels one by one as 2D+T; vertical aggregation will have to be done manually
for l in sorted(levels, reverse=True):
    t = r.readfield({'indicatorOfParameter':11,
                     'indicatorOfTypeOfLevel':100,
                     'level':l})
    temp_levels.append(t)
In [49]:
temp_levels
Out[49]:
[<epygram.fields.H2DField.H2DField at 0x7f333f53ba10>,
 <epygram.fields.H2DField.H2DField at 0x7f333f603c50>,
 <epygram.fields.H2DField.H2DField at 0x7f333f601b10>,
 <epygram.fields.H2DField.H2DField at 0x7f333f5f9d90>,
 <epygram.fields.H2DField.H2DField at 0x7f333f5f9610>,
 <epygram.fields.H2DField.H2DField at 0x7f333f4ce510>,
 <epygram.fields.H2DField.H2DField at 0x7f3373403810>,
 <epygram.fields.H2DField.H2DField at 0x7f333ebd80d0>,
 <epygram.fields.H2DField.H2DField at 0x7f333ebd8f50>,
 <epygram.fields.H2DField.H2DField at 0x7f333f601cd0>,
 <epygram.fields.H2DField.H2DField at 0x7f3373403750>,
 <epygram.fields.H2DField.H2DField at 0x7f333f601c50>,
 <epygram.fields.H2DField.H2DField at 0x7f333ebd8050>,
 <epygram.fields.H2DField.H2DField at 0x7f333ec85250>,
 <epygram.fields.H2DField.H2DField at 0x7f333ebd8dd0>]
In [50]:
temp_4D_data = numpy.concatenate([f.getdata(d4=True) for f in temp_levels], axis=1)
# axis = 1 because data in epygram fields is stored with order: (t,z,y,x)
In [51]:
temp_4D_data.shape
Out[51]:
(3, 15, 601, 801)
In [52]:
# then create a 4Dfield from a 2D+T one
temp_4d = temp_levels[0].deepcopy()  # now still 2D+T
# create an ad hoc vertical geometry, making a copy where we just change the list of levels
vgeom_4d = temp_4d.geometry.vcoordinate.footprint_clone(extra={'levels':sorted(levels, reverse=True)})
temp_4d.geometry.vcoordinate = vgeom_4d
In [53]:
# finally, set 4D data inside
temp_4d.setdata(temp_4D_data)
In [54]:
temp_4d.fid['netCDF'] = 'temperature'
out = epygram.formats.resource('out4D_manual.nc', 'w', fmt='netCDF')
out.writefield(temp_4d)
In [55]:
out.close()
out.open(openmode='r')
out.what()
### FORMAT: netCDF

NetCDF Global Attributes:
        made_with: u'epygram-1.2.10'
        Conventions: u'CF-1.6'
NetCDF dimension information:
        Name: time
                size: 3
        Name: Z
                size: 15
        Name: Y
                size: 601
        Name: X
                size: 801
NetCDF variable information:
        Name: time
                dimensions: (u'time',)
                size: 3
                type: dtype('float64')
                units: u'seconds since 2015-02-16 00:00:00'
        Name: Z
                dimensions: (u'Z',)
                size: 15
                type: dtype('float64')
                units: u'hPa'
                short_name: u'pressure'
        Name: longitude
                dimensions: (u'Y', u'X')
                size: 481401
                type: dtype('float64')
                units: u'degrees'
        Name: latitude
                dimensions: (u'Y', u'X')
                size: 481401
                type: dtype('float64')
                units: u'degrees'
        Name: temperature
                dimensions: (u'time', u'Z', u'Y', u'X')
                size: 21663045
                type: dtype('float64')
                units: u'K'