I want to consolidate channel 1, channel 31, the lats and the lons in a smaller file so that I don't have to store all 147 Mbytes of the original L1B file. In this notebook I'll write out a new h5 file following the h5py quick start tutorial: http://docs.h5py.org/en/latest/quick.html. To do this, I'll start from the satelliteI_h5 notebook and add a file writing section.
In [0]:
from __future__ import print_function
import os,site
import glob
import h5py
#
# add the lib folder to the path assuming it is on the same
# level as the notebooks folder
#
libdir=os.path.abspath('../lib')
site.addsitedir(libdir)
from modismeta_h5 import parseMeta
from h5dump import dumph5
In [0]:
ls
the glob function finds a file using a wildcard to save typing (google: python glob wildcard)
In [0]:
l1b_filename=glob.glob('../data/MYD02*.h5')[0]
print("found l1b file {}".format(l1b_filename))
geom_filename=glob.glob('../data/MYD03*.h5')[0]
print("found geom file {}".format(geom_filename))
In [0]:
l1b_file=h5py.File(l1b_filename)
geom_file=h5py.File(geom_filename)
Read the channel 31 radiance data from MODIS_SWATH_Type_L1B/Data Fields/EV_1KM_Emissive
In [0]:
print(l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['Band_1KM_Emissive'].shape)
print(l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['Band_1KM_Emissive'][...])
print(l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'].shape)
note that channel 31 occurs at index value 10
In [0]:
index31=10
the data is stored as unsigned, 2 byte integers which can hold values from 0 to $2^{16}$ - 1 = 65,535
In [0]:
chan31=l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'][index31,:,:]
print(chan31.shape,chan31.dtype)
In [0]:
chan31[:3,:3]
we need to apply a scale and offset to convert to radiance (the netcdf module did this for us automatically
$Data = (RawData - offset) \times scale$
this information is included in the attributes of each variable.
(see page 36 of the Modis users guide )
In [0]:
scale=l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'].attrs['radiance_scales'][index31]
offset=l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'].attrs['radiance_offsets'][index31]
In [0]:
chan31=(chan31 - offset)*scale
In [0]:
%matplotlib inline
histogram the calibrated radiances and show that they lie between 0-10 $W\,m^{-2}\,\mu m^{-1}\,sr^{-1}$
In [0]:
import matplotlib.pyplot as plt
out=plt.hist(chan31.flat)
Now do the same for Channel 1
From the Modis users guide (p. 25, table 3.3.2), we know that the 1 km version of channel 1 is called EV_250_Aggr1km_RefSB and is at index 0, so just get that channel, scale and offset for index 0
In [0]:
reflective=l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['EV_250_Aggr1km_RefSB'][0,:,:]
In [0]:
scale=l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['EV_250_Aggr1km_RefSB'].attrs['reflectance_scales']
offset=l1b_file['MODIS_SWATH_Type_L1B']['Data Fields']['EV_250_Aggr1km_RefSB'].attrs['reflectance_offsets']
chan1=(reflective - offset[0])*scale[0]
In [0]:
out=plt.hist(chan1.flat)
In [0]:
the_lon=geom_file['MODIS_Swath_Type_GEO']['Geolocation Fields']['Longitude'][...]
the_lat=geom_file['MODIS_Swath_Type_GEO']['Geolocation Fields']['Latitude'][...]
In [0]:
out_name="swath_output.h5"
f=h5py.File(out_name, "w")
dset = f.create_dataset("lattitude", the_lat.shape, dtype=the_lat.dtype)
dset[...]=the_lat[...]
dset = f.create_dataset("longitude", the_lon.shape, dtype=the_lon.dtype)
dset[...]=the_lon[...]
dset = f.create_dataset("channel1", chan1.shape, dtype=chan1.dtype)
dset[...]=chan1[...]
dset = f.create_dataset("channel31", chan31.shape, dtype=chan31.dtype)
dset[...]=chan31[...]
In [0]:
metadata=parseMeta(l1b_file)
metadata
Transfer all the metadata to the new hdf5 file as global attributes
In [0]:
for the_key in metadata.keys():
f.attrs[the_key]=metadata[the_key]
In [0]:
f.close()
Now reopen the output file and dump it to make sure it has all our stuff
In [0]:
f=h5py.File(out_name,'r')
dumph5(f)
f.close()