## "state33" data file obtained from: ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/daily/READMEimport gzipimport numpy as npimport matplotlib.pyplot as pltfid = gzip.open("state33.gz")month1 = 3month2 = 7M1,M2 = [],[]for line in fid:    stid = int(line[0:7])    if line[7:11]=="TMAX":        mo = int(line[17:19])        nday = int(line[20:22])        if mo not in [month1,month2]: continue        V,ii = [],26        for k in range(nday):            V.append(float(line[ii:ii+2]))            ii += 8        if mo==month1: M1.extend(V)        if mo==month2: M2.extend(V)M1 = np.array(np.sort(M1))M2 = np.array(np.sort(M2))M1 = M1[M1!=99]M2 = M2[M2!=99]M1 = (M1-M1.mean())/M1.std()M2 = (M2-M2.mean())/M2.std()Q1 = [M1[int(q*len(M1))] for q in np.arange(1,1000,dtype=np.float64)/1000]Q2 = [M2[int(q*len(M2))] for q in np.arange(1,1000,dtype=np.float64)/1000]P1 = [M1[int(q*len(M1))] for q in np.arange(1,10,dtype=np.float64)/10]P2 = [M2[int(q*len(M2))] for q in np.arange(1,10,dtype=np.float64)/10]plt.clf()plt.figure(figsize=(4,3.5))plt.axes([0.15,0.15,0.8,0.8])plt.grid(True)plt.plot(Q1, Q2, '-', color="gray", lw=3)plt.hold(True)plt.plot([-3,3], [-3,3], '-', color='black')plt.plot(P1, P2, 'o', mec='black', mfc='red')plt.xlim(-2,3)plt.ylim(-2,3)plt.xlabel("March", size=14)plt.ylabel("July", size=14)plt.savefig("ohio_temps_qq.pdf")plt.savefig("ohio_temps_qq.svg")