text: """
  NCL_stream_1.py
  ===============
  This script illustrates the following concepts:
     - Drawing a black-and-white streamline plot over a map
  
  See following URLs to see the reproduced NCL plot & script:
      - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/stream_1.ncl
      - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/stream_1_lg.png
  
  Note: The streamlines generated in this plot approximate, but do not exactly
  match those generated by the NCL version. This is primarily because the
  seed point generation algorithm used by NCL is opaque.
  """
  
  import cartopy.crs as ccrs
  import cartopy.feature as cfeature
  import geocat.datafiles as gdf
  import matplotlib.pyplot as plt
  ################################################################################
  # Import packages:
  import numpy as np
  import xarray as xr
  from geocat.viz import util as gvutil
  
  ################################################################################
  # Read in data:
  
  # Open a netCDF data file using xarray default engine and load the data into xarrays
  ds = xr.open_dataset(gdf.get('netcdf_files/uvt.nc'))
  # Extract a 2D horizontal slice from the first time step of the 3D U and V variables at the bottom level
  U = ds.U.isel(time=0, lev=0)
  V = ds.V.isel(time=0, lev=0)
  
  ################################################################################
  # Plot:
  
  # Generate figure (set its size (width, height) in inches)
  plt.figure(figsize=(16, 8))
  
  # Generate axes using Cartopy projection
  projection = ccrs.PlateCarree()
  ax = plt.axes(projection=projection)
  
  # Use global map
  ax.set_global()
  
  # Stream-plot the data
  # There is no Xarray streamplot function, yet. So need to call matplotlib.streamplot directly. Not sure why, but can't
  # pass xarray.DataArray objects directly: fetch NumPy arrays via 'data' attribute'
  ax.streamplot(U.lon.data,
                U.lat.data,
                U.data,
                V.data,
                linewidth=1,
                density=4,
                color='black',
                zorder=1)
  
  # Use geocat.viz.util convenience function to add minor and major tick lines
  gvutil.add_major_minor_ticks(ax, labelsize=16)
  
  # Use geocat.viz.util convenience function to make plots look like NCL plots by using latitude, longitude tick labels
  gvutil.add_lat_lon_ticklabels(ax)
  
  # Use geocat.viz.util convenience function to set axes tick values without calling two different matplotlib functions
  gvutil.set_axes_limits_and_ticks(ax,
                                   xticks=np.linspace(-180, 180, 13),
                                   yticks=np.linspace(-90, 90, 7))
  
  # Draw filled polygons for land
  ax.add_feature(cfeature.LAND, zorder=0, edgecolor='black', color='lightgray')
  
  # Use geocat.viz.util convenience function to add titles to left and right of the plot axis.
  gvutil.set_titles_and_labels(ax,
                               maintitle="Example of a streamline plot",
                               maintitlefontsize=22,
                               lefttitle=U.long_name,
                               lefttitlefontsize=18,
                               righttitle=U.units,
                               righttitlefontsize=18,
                               xlabel="",
                               ylabel="")
  
  # Show the plot
  plt.tight_layout()
  plt.show()
  