Merge pull request #569 from kafitzgerald/issue564

Update NCL_topo_1.py example to use the original dataset

Gallery/Topography/NCL_topo_1.py

@@ -0,0 +1,102 @@
+"""
+NCL_topo_1.py
+===============
+This script illustrates the following concepts:
+   - Drawing a topographic map using 5' data
+   - Drawing topographic data using GMT colormap
+   - Reading binary data
+   - Creating an xarray DataArray
+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
+"""
+
+###############################################################################
+# 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 geocat.datafiles as gdf
+import geocat.viz as gv
+
+###############################################################################
+# Read in data:
+
+# Open the dataset using numpy
+nlat = 2160
+nlon = 4320
+
+elevation_data = np.fromfile(gdf.get("binary_files/ETOPO5.DAT"),
+                             dtype='>i2').reshape((nlat, nlon))
+
+# Create numpy arrays for latitude and longitude
+lat = np.linspace(90, -90, nlat)
+lon = np.linspace(0, 360, nlon)
+
+# Create an xarray DataArray
+da = xr.DataArray(data=elevation_data,
+                  dims=["lat", "lon"],
+                  coords=dict(lat=(["lat"], lat, {
+                      "long_name": "latitude"
+                  }),
+                              lon=(["lon"], lon, {
+                                  "long_name": "longitude"
+                              })),
+                  name="elevation",
+                  attrs={"units": "m"})
+
+###############################################################################
+# 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 = da.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('ETOPO5', 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()