## "state33" data file obtained from: ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/daily/README

import gzip
import numpy as np
import matplotlib.pyplot as plt

fid = gzip.open("state33.gz")

month1 = 3
month2 = 7

M1,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")