In [1]:

    

from a301utils.a301_readfile import download
from a301lib.cloudsat import get_geo
import glob
import os
from pathlib import Path
import sys
import json
import numpy as np
import h5py
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap



rad_file='MYD021KM.A2006303.2220.006.2012078143305.h5'
geom_file='MYD03.A2006303.2220.006.2012078135515.h5'
lidar_file='2006303212128_02702_CS_2B-GEOPROF-LIDAR_GRANULE_P2_R04_E02.h5'
download(rad_file)
download(geom_file)
download(lidar_file)


    

    




MYD021KM.A2006303.2220.006.2012078143305.h5 already exists
and is 153363446 bytes
will not overwrite


MYD03.A2006303.2220.006.2012078135515.h5 already exists
and is 27751661 bytes
will not overwrite


2006303212128_02702_CS_2B-GEOPROF-LIDAR_GRANULE_P2_R04_E02.h5 already exists
and is 20991920 bytes
will not overwrite

In [2]:

    

lats,lons,date_times,prof_times,dem_elevation=get_geo(lidar_file)


    

In [3]:

    

from a301utils.modismeta_read import parseMeta
metadict=parseMeta(rad_file)
corner_keys = ['min_lon','max_lon','min_lat','max_lat']
min_lon,max_lon,min_lat,max_lat=[metadict[key] for key in corner_keys]


    

In [4]:

    

lon_hit=np.logical_and(lons>min_lon,lons<max_lon)
lat_hit = np.logical_and(lats>min_lat,lats< max_lat)
in_box=np.logical_and(lon_hit,lat_hit)
print("ground track has {} points, we've selected {}".format(len(lon_hit),np.sum(in_box)) )
box_lons,box_lats=lons[in_box],lats[in_box]


    

    




ground track has 37082 points, we've selected 2305

In [9]:

    

from a301lib.modis_reproject import make_projectname
reproject_name=make_projectname(rad_file)
reproject_path = Path(reproject_name)

if reproject_path.exists():
        print('using reprojected h5 file {}'.format(reproject_name))
else: #need to create reproject.h5 for channel 1
    channels='-c 1 4 3 31'
    template='python -m a301utils.modis_to_h5 {} {} {}'
    command=template.format(rad_file,geom_file,channels)
    if 'win' in sys.platform[:3]:
       print('platform is {}, need to run modis_to_h5.py in new environment'
          .format(sys.platform))
       print('open an msys terminal and run \n{}\n'.format(command))      
    else: #osx, so presample is available
        print('running \n{}\n'.format(command))
        out=os.system(command)
        the_size=reproject_path.stat().st_size
        print('generated reproject file for 4 channels, size is {} bytes'.format(the_size))


    

    




running 
python -m a301utils.modis_to_h5 MYD021KM.A2006303.2220.006.2012078143305.h5 MYD03.A2006303.2220.006.2012078135515.h5 -c 1 4 3 31

generated reproject file for 4 channels, size is 68407776 bytes

In [6]:

    

with h5py.File(reproject_name,'r') as h5_file:
    basemap_args=json.loads(h5_file.attrs['basemap_args'])
    chan1=h5_file['channels']['1'][...]
    geo_string = h5_file.attrs['geotiff_args']
    geotiff_args = json.loads(geo_string)

print('basemap_args: \n{}\n'.format(basemap_args))
print('geotiff_args: \n{}\n'.format(geotiff_args))


    

    




basemap_args: 
{'lat_0': 44.49704429255, 'urcrnrlat': 54.11341423457169, 'projection': 'laea', 'rsphere': [6378137.0, 6356752.314245179], 'llcrnrlon': -155.74610198251494, 'urcrnrlon': -117.16732627370857, 'llcrnrlat': 32.88085960245767, 'lon_0': -142.03761469569002}

geotiff_args: 
{'proj4_string': '+units=m +proj=laea +datum=WGS84 +lon_0=-142.04 +lat_0=44.50', 'proj_id': 'laea', 'adfgeotransform': [-1282221.350437168, 1300.1525690949175, 0, 1324541.1860649737, 0, -1300.6679472554665], 'width': 2215, 'height': 1930}

In [7]:

    

%matplotlib inline
from matplotlib import cm
from matplotlib.colors import Normalize

cmap=cm.autumn  #see http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
cmap.set_over('w')
cmap.set_under('b',alpha=0.2)
cmap.set_bad('0.75') #75% grey

plt.close('all')
fig,ax = plt.subplots(1,1,figsize=(14,14))
#
# set up the Basemap object
#
basemap_args['ax']=ax
basemap_args['resolution']='c'
bmap = Basemap(**basemap_args)
#
# transform the ground track lons/lats to x/y
#
cloudsatx,cloudsaty=bmap(box_lons,box_lats)
#
# plot as blue circles
#
bmap.plot(cloudsatx,cloudsaty,'bo')
#
# now plot channel 1
#
num_meridians=180
num_parallels = 90
col = bmap.imshow(chan1, origin='upper',cmap=cmap, vmin=0, vmax=0.4)

lon_sep, lat_sep = 5,5
parallels = np.arange(-90, 90, lat_sep)
meridians = np.arange(0, 360, lon_sep)
bmap.drawparallels(parallels, labels=[1, 0, 0, 0],
                       fontsize=10, latmax=90)
bmap.drawmeridians(meridians, labels=[0, 0, 0, 1],
                       fontsize=10, latmax=90)
bmap.drawcoastlines()
colorbar=fig.colorbar(col, shrink=0.5, pad=0.05,extend='both')
colorbar.set_label('channel1 reflectivity',rotation=-90,verticalalignment='bottom')
_=ax.set(title='vancouver')


    

In [8]:

    

groundtrack_name = reproject_name.replace('reproject','groundtrack')
print('writing groundtrack to {}'.format(groundtrack_name))
box_times=date_times[in_box]
#
# h5 files can't store dates, but they can store floating point
# seconds since 1970, which is called POSIX timestamp
#
timestamps = [item.timestamp() for item in box_times]
timestamps= np.array(timestamps)
with h5py.File(groundtrack_name,'w') as groundfile:
    groundfile.attrs['cloudsat_filename']=lidar_file
    groundfile.attrs['modis_filename']=rad_file
    groundfile.attrs['reproject_filename']=reproject_name
    dset=groundfile.create_dataset('cloudsat_lons',box_lons.shape,box_lons.dtype)
    dset[...] = box_lons[...]
    dset.attrs['long_name']='cloudsat longitude'
    dset.attrs['units']='degrees East'
    dset=groundfile.create_dataset('cloudsat_lats',box_lats.shape,box_lats.dtype)
    dset[...] = box_lats[...]
    dset.attrs['long_name']='cloudsat latitude'
    dset.attrs['units']='degrees North'
    dset= groundfile.create_dataset('cloudsat_times',timestamps.shape,timestamps.dtype)
    dset[...] = timestamps[...]
    dset.attrs['long_name']='cloudsat UTC datetime timestamp'
    dset.attrs['units']='seconds since Jan. 1, 1970'


    

    




writing groundtrack to MYD021KM.A2006303.2220.006.groundtrack.h5

In [ ]: