Merge pull request #453 from NCAR/topo_1

Topo 1

Gallery/Topography/NCL_topo_1.py

@@ -0,0 +1,96 @@
+"""
+NCL_topo_1.py
+===============
+This script illustrates the following concepts:
+   - Drawing a topographic map using 1' data
+   - Drawing topographic data using GMT colormap
+   - Using wget to download a dataset
+
+See following URLs to see the reproduced NCL plot & script:
+    - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/topo_1.ncl
+    - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/topo_1_lg.png
+
+Note:
+    In the original NCL script, the ETOPO5 dataset was used. For this example,
+    we use the most current version of this data, ETOPO1.
+"""
+
+###############################################################################
+# Import packages:
+
+import matplotlib.pyplot as plt
+import xarray as xr
+import numpy as np
+import cartopy.crs as ccrs
+from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter
+import cmaps
+import requests
+
+import geocat.viz as gv
+
+###############################################################################
+# Read in data:
+
+# Use requests module to retrieve and download NOAA elevation data
+url = 'https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/ETOPO1_Ice_g_gmt4.grd.gz'
+req = requests.get(url)
+filename = 'ETOPO1_Ice_g_gmt4.grd.gz'
+with open(filename, 'wb') as output_file:
+    output_file.write(req.content)
+
+# Open the dataset using xarray
+ds = xr.open_dataset(filename)
+
+# Select the elevation data
+ds = ds.z
+
+###############################################################################
+# Plot
+
+# Generate figure and set size
+plt.figure(figsize=(12, 10))
+
+# Generate axes, using Cartopy
+projection = ccrs.PlateCarree()
+ax = plt.axes(projection=projection)
+
+# Add coastlines
+ax.coastlines(zorder=10)
+
+# Plot the elevation data
+elev = ds.plot.imshow(ax=ax,
+                      transform=projection,
+                      cmap=cmaps.GMT_relief,
+                      add_colorbar=False)
+
+# Add colorbar
+cbar = plt.colorbar(ax=ax, mappable=elev, orientation='horizontal', pad=0.1)
+cbar.ax.tick_params(
+    size=0,
+    labelsize=14)  # Remove the tick marks from the colorbar, set label size
+cbar.ax.xaxis.set_tick_params(pad=10)
+
+# Use geocat-viz utility function to format major and minor tick marks
+gv.add_major_minor_ticks(ax, labelsize=14)
+
+# Add title, remove default x and y labels
+plt.title('ETOPO1', fontweight='bold', fontsize=23, pad=15)
+plt.xlabel("")
+plt.ylabel("")
+
+# Use geocat-viz utility function to format x and y tick labels
+gv.set_axes_limits_and_ticks(ax,
+                             xlim=[-180, 180],
+                             ylim=[-90, 90],
+                             xticks=np.arange(-180, 181, 30),
+                             yticks=np.arange(-90, 91, 30))
+
+# Use geocat-viz utility function to add lat/lon formatting for tick labels
+gv.add_lat_lon_ticklabels(ax)
+
+# Remove degree symbol from lat/lon labels
+ax.yaxis.set_major_formatter(LatitudeFormatter(degree_symbol=''))
+ax.xaxis.set_major_formatter(LongitudeFormatter(degree_symbol=''))
+
+# Show the plot
+plt.show()