<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Ploting points and lines on a reproject raster
Here's how to draw a line between seattle and portland given the raster modis image of Vancouver
produced by modis_to_h5.py
Note that this assumes channel 1 is named '1', not 'chan1' as in the older version of modis_to_h5.py
1. read in the chan1 read in the basemap argument string and turn it  into a dictionary of basemap arguments using json.loads
Step1:  2. Plot chan1 as in resample2.ipynb, and add vancouver and portland points with a line between them
Step2:  What is the distance between Vancouver and Portland?
pyproj.Geod defines a geoid using the major and minor axes from our Basemap datum and calculates
the azimuthal angles and distance between two points along a great circle rounte (i.e. shortest distance along the surface of the WGS84 geoid)
<ASSISTANT_TASK:>
Python Code:
import h5py
from a301utils.a301_readfile import download
from mpl_toolkits.basemap import Basemap
from matplotlib import pyplot as plt
import json
import numpy as np
rad_file=' MYD021KM.A2016217.1915.006.2016218155919.h5'
geom_file='MYD03.A2016217.1915.006.2016218154759.h5'
download(rad_file)
data_name='MYD021KM.A2016224.2100.006_new.reproject.h5'
download(data_name)
with h5py.File(data_name,'r') as h5_file:
    basemap_args=json.loads(h5_file.attrs['basemap_args'])
    chan1=h5_file['channels']['1'][...]
print(basemap_args)
Explanation: Ploting points and lines on a reproject raster
Here's how to draw a line between seattle and portland given the raster modis image of Vancouver
produced by modis_to_h5.py
Note that this assumes channel 1 is named '1', not 'chan1' as in the older version of modis_to_h5.py
1. read in the chan1 read in the basemap argument string and turn it  into a dictionary of basemap arguments using json.loads
End of explanation
%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)
num_meridians=180
num_parallels = 90
vmin=None; vmax=None
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')
#
# now use the basemap object to project the portland and vancouver
# lon/lat coords into the xy lambert coordinate system
#
#  remember what the asterisk * argument expansion does:
#  if I have a list A=[a,b] then fun(*A) is the same as fun(a,b)
#
#
vancouver_lon_lat=[-123.1207,49.2827]
portland_lon_lat=[-122.6765,45.5231]
#
# get the xy coords
#
van_xy = bmap(*vancouver_lon_lat)
portland_xy = bmap(*portland_lon_lat)
#
#  draw a blue circle for van and
#  a green circle for portland
#
bmap.plot(*van_xy,'bo',markersize=15)
bmap.plot(*portland_xy,'go',markersize=15)
#
# connect them with a cyan line
#
xcoords=[van_xy[0],portland_xy[0]]
ycoords=[van_xy[1],portland_xy[1]]
_ = bmap.plot(xcoords,ycoords,'c-',linewidth=5)
Explanation: 2. Plot chan1 as in resample2.ipynb, and add vancouver and portland points with a line between them
End of explanation
import pyproj
great_circle=pyproj.Geod(a=bmap.rmajor,b=bmap.rminor)
azi12,azi21,distance=great_circle.inv(vancouver_lon_lat[0],vancouver_lon_lat[1],
                                                        portland_lon_lat[0],portland_lon_lat[1])
print('Vancouver to Portland -- great circle is: {:5.2f} km'.format(distance/1.e3))
Explanation: What is the distance between Vancouver and Portland?
pyproj.Geod defines a geoid using the major and minor axes from our Basemap datum and calculates
the azimuthal angles and distance between two points along a great circle rounte (i.e. shortest distance along the surface of the WGS84 geoid)
End of explanation
<END_TASK>