5.4. Make 4D fields (to netCDF)

5.4.1. with integrated functionalities: meta-resources

[1]:
import os
import epygram
epygram.init_env()
workdir = epygram.config.userlocaldir + '/notebooks_data'
os.chdir(workdir)
import numpy
[2]:
files = ['advanced_examples/ICMSHAROM+{:0>4}'.format(i) for i in range(0,2+1)]
[3]:
# 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
[4]:
print(t.data.shape)
t.sp2gp()
print(t.data.shape)
(3, 60, 423256)
(3, 60, 720, 750)
[5]:
rh = r.readfield({'discipline':0, 'parameterCategory':1, 'parameterNumber':1,
                  'typeOfFirstFixedSurface':103, 'level':2})  # relative humidity at 2m, as GRIB2
[6]:
out = epygram.formats.resource('out4D.nc', 'w', fmt='netCDF')  # open the output netCDF
[7]:
t.fid['netCDF'] = 'temperature'
out.writefield(t)
[8]:
rh.fid['netCDF'] = 'relative_humidity'
out.writefield(rh)
# [2018/12/19-13:45:08][epygram.formats.netCDF][writefield:1100][INFO]: assume lons/lats match.
# [2018/12/19-13:45:08][epygram.formats.netCDF][writefield:1205][INFO]: assume projection parameters match.
[9]:
out.close()
out.open(openmode='r')
out.what()
### FORMAT: netCDF

NetCDF Global Attributes:
        made_with: u'epygram-1.3.5'
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')
                _FillValue: -999999.90000000002
                grid_mapping: u'Projection_parameters'
                vertical_grid: u'Z'
        Name: relative_humidity
                dimensions: (u'time', u'Y', u'X')
                size: 1620000
                type: dtype('float64')
                _FillValue: -999999.90000000002
                grid_mapping: u'Projection_parameters'

5.4.2. or do it manually

[10]:
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']
[11]:
# try with Combine-Levels meta-resource
r = epygram.resources.meta_resource(files[0], 'r', 'CL')
[12]:
t = r.readfield({'indicatorOfParameter':11,
                 'indicatorOfTypeOfLevel':100})  # temperature on Pressure levels, as GRIB1

AssertionErrorTraceback (most recent call last)
<ipython-input-12-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)
    160         :param getdata: if False, do not read data but only metadata
    161         """
--> 162         result = self.readfields(handgrip=handgrip, getdata=getdata)
    163         if len(result) != 1:
    164             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)
    173         """
    174         fieldset = FieldSet()
--> 175         cont = self._create_list()
    176         for fid in self.listfields(select=handgrip):
    177             fid = FPDict(fid)  # footprints purpose

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

AssertionError: Not able to combine levels if fields do not have 'generic' fids
[ ]:
# This did not worked because 'CL' meta-resource need generic/GRIB2 => let's do it manually
[13]:
r = epygram.resources.meta_resource(files, 'r', 'MV')  # this works anyway
[14]:
t = r.readfield({'indicatorOfParameter':11,
                 'indicatorOfTypeOfLevel':100,
                 'level':850})  # temperature, as GRIB1
[15]:
t.data.shape
[15]:
(3, 601, 801)
[16]:
# 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)
[17]:
# 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]
[18]:
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)
[19]:
temp_levels
[19]:
[<epygram.fields.H2DField.H2DField at 0x7f5d780233d0>,
 <epygram.fields.H2DField.H2DField at 0x7f5d6572c7d0>,
 <epygram.fields.H2DField.H2DField at 0x7f5d6572cd90>,
 <epygram.fields.H2DField.H2DField at 0x7f5d6579df50>,
 <epygram.fields.H2DField.H2DField at 0x7f5d90777990>,
 <epygram.fields.H2DField.H2DField at 0x7f5d6579ded0>,
 <epygram.fields.H2DField.H2DField at 0x7f5d6572c210>,
 <epygram.fields.H2DField.H2DField at 0x7f5d6583d290>,
 <epygram.fields.H2DField.H2DField at 0x7f5d6583d690>,
 <epygram.fields.H2DField.H2DField at 0x7f5d657b1290>,
 <epygram.fields.H2DField.H2DField at 0x7f5d6583d210>,
 <epygram.fields.H2DField.H2DField at 0x7f5d657b1a10>,
 <epygram.fields.H2DField.H2DField at 0x7f5d65848c50>,
 <epygram.fields.H2DField.H2DField at 0x7f5d65848e10>,
 <epygram.fields.H2DField.H2DField at 0x7f5d657b1650>]
[20]:
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)
[21]:
temp_4D_data.shape
[21]:
(3, 15, 601, 801)
[22]:
# 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
[23]:
# finally, set 4D data inside
temp_4d.setdata(temp_4D_data)
[24]:
temp_4d.fid['netCDF'] = 'temperature'
out = epygram.formats.resource('out4D_manual.nc', 'w', fmt='netCDF')
out.writefield(temp_4d)
[25]:
out.close()
out.open(openmode='r')
out.what()
### FORMAT: netCDF

NetCDF Global Attributes:
        made_with: u'epygram-1.3.5'
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')
                _FillValue: -999999.90000000002
                units: u'K'
[ ]: