Dependencies | Version |
---|---|
SatPy | 0.9.0 |
This example loads multiple VIIRS SDR granules of band I04 data, resamples them to a known AreaDefinition
, and uses matplotlib to plot the data on a map. Although cartopy is not imported directly, we use pyresample's to_cartopy_crs()
method to create a cartopy CRS object and use that in the plotting.
This example should act as a starting point for plotting data processed with SatPy. Additional information can be found on XArray's plotting documentation as well as Cartopy's documentation for help building more complex plots.
Using the glob
module makes it easy to get a series of data files. After creating the Scene object, we load the I04 band data.
from satpy import Scene
from glob import glob
import matplotlib.pyplot as plt
filenames = glob('/data/data/viirs/conus_day/*t180*.h5')
scn = Scene(reader='viirs_sdr', filenames=filenames)
scn.load(['I04'])
Next we put the VIIRS data on a high resolution uniform grid. We could have also used any other area definition. For polar-orbiting satellite data it is best to resample the data first before putting it on a plot due to the size of the lon/lat arrays that would need to be used. For geo-stationary satellite data, which is typically already gridded, the source data's area can be used directly in the plotting steps below.
my_area = scn['I04'].attrs['area'].compute_optimal_bb_area({'proj': 'lcc', 'lon_0': -95., 'lat_0': 25., 'lat_1': 25., 'lat_2': 25.})
new_scn = scn.resample(my_area)
/Users/davidh/anaconda/envs/polar2grid_py36/lib/python3.6/site-packages/dask/local.py:271: RuntimeWarning: invalid value encountered in greater_equal return func(*args2) /Users/davidh/anaconda/envs/polar2grid_py36/lib/python3.6/site-packages/dask/local.py:271: RuntimeWarning: invalid value encountered in less_equal return func(*args2)
To plot the data using cartopy we first convert the pyresample AreaDefinition to a cartopy CRS and then use that to create matplotlib axes.
crs = new_scn['I04'].attrs['area'].to_cartopy_crs()
ax = plt.axes(projection=crs)
ax.coastlines()
ax.gridlines()
ax.set_global()
plt.imshow(new_scn['I04'], transform=crs, extent=crs.bounds, origin='upper')
cbar = plt.colorbar()
cbar.set_label("Kelvin")
plt.show()
/Users/davidh/repos/git/pyresample/pyresample/_cartopy.py:36: UserWarning: 'cartopy' >= 0.17 required for better 'from_proj' functionality. warnings.warn("'cartopy' >= 0.17 required for better 'from_proj' "
![title](Cartopy Plot 1.png)