# IPython session run with --pylab=inline; "import matplotlib.pyplot as plt" was run during startup
import numpy as np
import pandas as pd

# ERDDAP data request URL (decided in advance), selecting only a subset of two variables
# Returns a 'csvp' table with a single header row that concatenates variable names and units
pss2tags_url = 'http://MYSERVER/erddap/tabledap/otnnepPSS2AnTags.csvp?length,weight'

# read animal tag data into pandas dataframe, then rename columns to more 'robust' names (no spaces or parentheses)
resppd = pd.read_csv(pss2tags_url)
resppd.rename(columns={'length (meters)': 'length_m', 'weight (kg)': 'weight_kg'}, inplace=True)
resppd

# plot animal weight vs length, after converting units
x=resppd.length_m * 100
y=resppd.weight_kg * 1000
fig1 = plt.figure()
plt.plot(x,y,'.r')
plt.xlabel('length (cm)')
plt.ylabel('weight (gram)')

# Process data to generate 2D histogram
nbins = 60
H, xedges, yedges = np.histogram2d(x,y,bins=nbins)
# H needs to be rotated and flipped
H = np.rot90(H)
H = np.flipud(H)
# Mask pixels with a value of zero
Hmasked = np.ma.masked_where(H==0,H)

# Plot 2D histogram using pcolor
# This plot is similar to the previous one (eg, same x & y axes), but much more informative
fig2 = plt.figure()
plt.pcolormesh(xedges,yedges,Hmasked)
plt.xlabel('length (cm)')
plt.ylabel('weight (gram)')
cbar = plt.colorbar()
cbar.ax.set_ylabel('Counts')