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



In [1]:

    
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)









    



MYD021KM.A2016224.2100.006_new.reproject.h5 already exists, need to either delete or rename for new download
{'urcrnrlat': 55.64793157078861, 'rsphere': [6378137.0, 6356752.314245179], 'lat_0': 46.12922265815655, 'llcrnrlon': -136.53911711598428, 'llcrnrlat': 34.43216712325068, 'projection': 'laea', 'lon_0': -122.60317657850351, 'urcrnrlon': -96.49588633620655}

2. Plot chan1 as in resample2.ipynb, and add vancouver and portland points with a line between them



In [2]:

    
%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)

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)



In [5]:

    
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))









    



Vancouver to Portland -- great circle is: 419.33 km



In [ ]: