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'
[ ]: