import numpy as np
import matplotlib.cbook as cbook
import matplotlib.colors as colors
import matplotlib.pyplot as plt
# Load the example data
dem = cbook.get_sample_data('topobathy.npz', np_load=True)
# Compute elevations in meters instead of feet
topo = dem['topo']*0.3048
longitude = dem['longitude']
latitude = dem['latitude']
colors_undersea = plt.cm.terrain(np.linspace(0, 0.17, 256))
colors_land = plt.cm.terrain(np.linspace(0.25, 1, 256))
all_colors = np.vstack((colors_undersea, colors_land))
terrain_map = colors.LinearSegmentedColormap.from_list('terrain_map', all_colors)
divnorm = colors.TwoSlopeNorm(vmin=-200., vcenter=0, vmax=1200.)
# Plot our map and tell Matplotlib to use our new normalized colormap
pcm = plt.pcolormesh(longitude, latitude, topo, rasterized=True, norm=divnorm, cmap=terrain_map, shading='auto')
# Create a colorbar
cb = plt.colorbar(pcm, shrink=0.8)
cb.set_ticks([-200, 0, 400, 800, 1200])
cb.set_label('Elevation [m]')
plt.gca().set_aspect(1 / np.cos(np.deg2rad(49)))
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.tight_layout()
plt.savefig('colormap.png')