###################################################################################################
# There are multipele codes in here. Codes are not in order
###################################################################################################

############################################
# Read in crop yield data
# Finds lat and lon of counties
############################################
from pylab import *
import csv
from math import sqrt
from sys import exit
import pickle
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from mpl_toolkits.basemap import Basemap
import os
from scipy.stats import norm
import matplotlib as mpl
from matplotlib.patches import Polygon
import random

# cd Documents/Science_Fair_2017_Crop_Yields/

# the number of everything
nCounties=3143
nCrop=3
nYears=116
nStates=57

#files to open
f=open('data/crop_yields.csv','r')
ft=open('data/crop_totals.csv','r')
floc=open('data/countylatlon.csv')

presentGrowingCounties=pickle.load(open('pickle_files/presentGrowingCounties.p','rb'))

# initialize variables
iCounty=0 # a counter of the number of counties
cIndex=-9999*ones(shape=(nStates,850))
cropYield=-9999.*ones(shape=(nCounties,nYears,nCrop))
cropTotal=-9999.*ones(shape=(nCounties,nYears,nCrop))
countyName=[]
lat=-9999.*ones(shape=(nCounties))
lon=-9999.*ones(shape=(nCounties))
cIndexState={} # give it cIndex, returns state
highTotalYield=zeros(shape=(nCounties,nCrop),dtype=bool)

#                              #
# reads in county lat and lons #
#                              #
for line in floc:
    tmp=line.split(',')
    
    stateID=float(tmp[1])
    countyID=float(tmp[3])
    countyName.append(tmp[2])
    
    s=stateID
    c=countyID
        
    cIndex[s,c]=iCounty
    cIndexState[cIndex[s,c]]=s
    
    lat[cIndex[s,c]]=tmp[5]
    lon[cIndex[s,c]]=tmp[6]
    iCounty+=1
floc.close()
j=-1
countyName=[''  for x in xrange(3143)]

#                          #
# reads in crop yield data #
#                          #
for line in f:
    j+=1
    line=line.translate(None, ',')
    tmp=line.split('""')
    
    if tmp[9]=='OTHER (COMBINED) COUNTIES':
        continue
    
    year=float(tmp[1])
    stateID=float(tmp[6])
    countyID=float(tmp[10])
    crop=tmp[15]
    yd=float(tmp[19])
    
    if crop=='CORN':
        cropID=0
    if crop=='SOYBEANS':
        cropID=1
    if crop=='RICE':
        cropID=2
    
    y=year-1900
    s=stateID
    c=countyID
    cp=cropID
    
    if cIndex[s,c]==-9999:
        continue
    cropYield[cIndex[s,c],y,cp]=yd 
    countyName[int(cIndex[s,c])]=tmp[9]

f.close()

j=0
#                           #
# reads in total yield data #
#                           #
for line in ft:
    j+=1
    line=line.translate(None, ',')
    tmp=line.split('""')
    
    if tmp[9]=='OTHER (COMBINED) COUNTIES':
        continue
    
    year=float(tmp[1])
    stateID=float(tmp[6])
    countyName.append(tmp[9])
    countyID=float(tmp[10])
    crop=tmp[15]
    yd=float(tmp[19])
    
    if crop=='CORN':
        cp=0
    if crop=='SOYBEANS':
        cp=1
    if crop=='RICE':
        cp=2
        
    y=year-1900
    s=stateID
    c=countyID
    
    if cIndex[s,c]==-9999:
        continue
    cropTotal[cIndex[s,c],y,cp]=yd 
ft.close()

for cp in range(nCrop):
    cutoff=.1*(np.amax(cropTotal[:,115,cp]))
    for icity in range(3143):
        total=average(cropTotal[icity,105:116,cp])
        if total>=cutoff:
            highTotalYield[icity,cp]=True     
            
#pickle.dump(highTotalYield,open('pickle_files/highTotalYield.p','wb'))  
        
#pickle.dump(cropTotal,open('pickle_files/cropTotal.p','wb'))

# create the map
map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
        projection='lcc',lat_1=33,lat_2=45,lon_0=-95)

# load the shapefile, use the name 'states'
map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)

ax = plt.gca() # get current axes instance

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################

# collect the state names from the shapefile attributes so we can
# look up the shape obect for a state by it's name
j=0 #counter for how many times loop goes through
g=0 
cmapArray=plt.cm.jet(arange(256))
cmin=0
y1=0
y2=256
year=115
cp=2
cmax=.6*(np.amax(cropTotal[:,115,cp]))
cmin=np.amin(cropTotal[:,115,cp])
#cmax=np.amax(cropYield[:,115,cp])
#cmax=8

for shape_dict in map.states_info:
    seg = map.states[j]
    
    s=int(shape_dict['STATEFP'])
    c=int(shape_dict['COUNTYFP'])
    
    if s==72 or cIndex[s,c]==-9999 or presentGrowingCounties[cIndex[s,c],cp]==0:
        j+=1
        continue

    Yield=average(cropTotal[cIndex[s,c],105:116,cp])
    
    if Yield<=0:
        j+=1
        continue
        
    x=Yield
    y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
    icmap=min(255,int(round(y,1)))
       
    poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
    ax.add_patch(poly)
       
    j+=1
    g+=1
plt.show() 
    
# pickle the values
pickle.dump(lat,open('pickle_files/county_lats.p','wb'))
pickle.dump(lon,open('pickle_files/county_lons.p','wb'))
pickle.dump(cropYield,open('pickle_files/cropYield.p','wb'))
pickle.dump(cropTotal,open('pickle_files/cropTotal.p','wb'))
pickle.dump(cIndex,open('pickle_files/cIndex.p','wb'))
pickle.dump(cIndexState,open('pickle_files/cIndexState.p','wb'))  
pickle.dump(countyName,open('pickle_files/countyName.p','wb'))from pylab import *
import csv
from math import sqrt
from sys import exit
import numpy as np
import pickle
import os

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################

countyLat=pickle.load(open('pickle_files/county_lats.p','rb'))
countyLon=pickle.load(open('pickle_files/county_lons.p','rb'))
cmiplat=pickle.load(open('pickle_files/cmiplatall.p','rb'))
cmiplon=pickle.load(open('pickle_files/cmiplonall.p','rb'))
exit()
# initialize variables
nearestStationLat=-9999*ones(shape=(3143))
nearestStationLon=-9999*ones(shape=(3143))
latIndex=-9999*ones(shape=(3143))
lonIndex=-9999*ones(shape=(3143))

f=open('written_files/model_nearest_stations.txt','w')
zero=0

for c in range(3143):
    print c, ' of 3143'
    lat=countyLat[c]
    lon=countyLon[c]
    shortestDist=.15
    shortestDist2=.15
    b=0
    g=0
    for ilat in range(len(cmiplat)):
        modelLat=cmiplat[ilat]
        for ilon in range(len(cmiplon)):
            modelLon=cmiplon[ilon]-360
            
            # find distance from county center to weather station
            dist=abs(sqrt((lon-modelLon)**2+(lat-modelLat)**2))
                
            
            if dist>.1: # only look at the close ones
                b+=1
                continue
            
            if dist<shortestDist: # if its closer than the last closest one
                g+=1
                # make this one the closest one
                shortestDist=dist
                nearestStationLat[c]=modelLat
                nearestStationLon[c]=modelLon
                latIndex[c]=ilat
                lonIndex[c]=ilon
                         
    if g==0:
        zero+=1
        shortestDist=-9999
        nearestStationLat[c]=-9999
        nearestStationLon[c]=-9999
        latIndex[c]=-9999
        lonIndex[c]=-9999
        
    # write the closest station into a csv file      
    f.write(str(nearestStationLat[c])+','+str(nearestStationLon[c])+','+
        str(int(latIndex[c]))+','+str(int(lonIndex[c]))+'\n')

        
f.close()
            
###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################            
            
            
################################################################
# Finds two closest weather statioins to center of each county
################################################################
from pylab import *
import csv
from math import sqrt
from sys import exit
import pickle
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from mpl_toolkits.basemap import Basemap
import os
from scipy.stats import norm
import matplotlib as mpl
from matplotlib.patches import Polygon
import random

# cd Documents/Science_Fair_2017_Crop_Yields/

## Number of Years Wanted ##
years=45

###############################################
# Read in Station Names
###############################################
filename='data/noaa_daily_data/ghcnd-stations.txt'
openStations=open(filename,'r')

cityName={}  #initialize a dictionary

for line in openStations:
    tmp=line[0:85]
    station=tmp[0:11]
    cityName[station]=tmp[41:70]
       
###############################################
# Read in Start and End Date Data
###############################################

# load the pickle files
countyLat=pickle.load(open('pickle_files/county_lats.p','rb'))
countyLon=pickle.load(open('pickle_files/county_lons.p','rb'))

#the inventory file to read in
fileOpen=open('data/noaa_daily_data/ghcnd-inventory.txt','r') 

# a file I am going to write with Tmax for cities over however many years
f=open('written_files/tmax_cities_over_'+str(years)+'_years.txt','w')

# a file I am going to write with precip for cities over however many years
fprecip=open('written_files/precip_cities_over_'+str(years)+'_years_no.txt','w')

# a file I am going to write with the 2 closest stations to the county centers
fnearestStation=open('written_files/tmax_nearest_stations.txt','w')

fnearestStation2=open('written_files/tmax_second_nearest_stations.txt','w')

#fnearestStationsPrecip=open('written_files/precip_nearest_stations.txt','w')


#                      #
# initialize variables #
#                      #
k=0 #a counter of how many times the loop goes through
maxk=571434 # the maximum number of stations
startYear=zeros(shape=(maxk))
endYear=zeros(shape=(maxk))
length=zeros(shape=(maxk))
var=zeros(shape=(maxk))
stationLat=zeros(shape=(maxk))
stationLon=zeros(shape=(maxk))
stationLatTmax=zeros(shape=(50000))
stationLonTmax=zeros(shape=(50000))
kPlot=0
kPrecipPlot=0
stationLatPrecip=zeros(shape=(50000))
stationLonPrecip=zeros(shape=(50000))
longCity=zeros(shape=(maxk),dtype=bool)
stationCode=[]

#                                                                              #
# read in each the lat, lon, station name, startYear, endYear for each station #
#                                                                              #
for line in fileOpen: # loop through every station
    tmp=line[0:45]
    station=tmp[0:11]
    stationLat[k]=tmp[11:20]
    stationLon[k]=tmp[21:30]
    variable=tmp[31:35]
    startYear[k]=float(tmp[36:40])
    endYear[k]=float(tmp[41:45])
    
    stationCode.append(station)
    
    if variable=='TMAX':
        var[k]=0
    elif variable=='TMIN':
        var[k]=1
    else:
        var[k]=-1
    
    
    length[k]=endYear[k]-startYear[k] # find how many years of data the station has
            
    # find the cities that have atleast years worth of data, the var is Tmax, and ends after 2010
    if length[k]>=years and var[k]==0 and endYear[k]>2010: 
        longCity[k]=True
        stationLatTmax[kPlot]=stationLat[k] # put the latitude of the city into an array
        stationLonTmax[kPlot]=stationLon[k] # put the longitude of the city into an array
        kPlot+=1 # a counter of how many cities
        # write the station, lat, lon, var, startYear, endYear, and place where located into a document
        f.write(tmp+' '+cityName[station]+'\n')
     
    # find the cities that have atleast years worth of data, the var is precip, and ends after 2010    
    if length[k]>=years and var[k]==2 and endYear[k]>2010:
        # write the station, lat, lon, var, startYear, endYear, and place where located into a document
        fprecip.write(tmp+' '+cityName[station]+'\n')
        stationLatPrecip[kPrecipPlot]=stationLat[k] # put the latitude of the city into an array
        stationLonPrecip[kPrecipPlot]=stationLon[k] # put the longitude of the city into an array
        kPrecipPlot+=1 # a counter of how many cities
               
    k+=1 # add one to the counter of how many times the loop has gone through

f.close()
fprecip.close()

#print kPlot # print number of cities written into the Tmax file
#print kPrecipPlot # print number of cities written into the precip file

# initialize variables
nearestStationLat=-9999*ones(shape=(3143))
nearestStationLon=-9999*ones(shape=(3143))
nearestStationIndex=[-9999]*3143

nearestStationLat2=-9999*ones(shape=(3143))
nearestStationLon2=-9999*ones(shape=(3143))
nearestStationIndex2=[-9999]*3143

for c in range(3143):
    print c, ' of 3143'
    lat=countyLat[c]
    lon=countyLon[c]
    for w in range(maxk):
        if w==0: # only if it is the first time running through
            shortestDist=dist
        if lat==-9999:
            continue
        if length[w]<years or endYear[w]<2010: # dont want weather stations too short
            continue   
        if var[w]!=0: #only want tmax
            continue 
        wsLat=stationLat[w]
        wsLon=stationLon[w]
        
        # find distance from county center to weather station
        dist=abs(sqrt((lon-wsLon)**2+(lat-wsLat)**2))
        
        if dist>1: # only look at the close ones
            continue
        
        if dist<shortestDist: # if its closer than the last closest one
            # make the last closest one the second to closest one
            shortestDist2=shortestDist
            nearestStationLat2[c]=nearestStationLat[c]
            nearestStationLon2[c]=nearestStationLon[c]
            nearestStationIndex2[c]=nearestStationIndex[c]
            
            # make this one the closest one
            shortestDist=dist
            nearestStationLat[c]=wsLat
            nearestStationLon[c]=wsLon  
            nearestStationIndex[c]=w
              
    # write the closest and second closest stations into a csv file      
    fnearestStation.write(str(stationCode[nearestStationIndex[c]])+','+
        '1,'+
        str(nearestStationLat[c])+','+
        str(nearestStationLon[c])+','+
        str(var[nearestStationIndex[c]])+','+
        str(startYear[nearestStationIndex[c]])+','+
        str(endYear[nearestStationIndex[c]])+','+
        cityName[stationCode[nearestStationIndex[c]]]+','+
        str(stationCode[nearestStationIndex2[c]])+','+
        '2,'+
        str(nearestStationLat2[c])+','+
        str(nearestStationLon2[c])+','+
        str(var[nearestStationIndex2[c]])+','+
        str(startYear[nearestStationIndex2[c]])+','+
        str(endYear[nearestStationIndex2[c]])+','+
        cityName[stationCode[nearestStationIndex2[c]]] +'\n')

fnearestStation.close()        
  
#plot the closest stations
figure(1,figsize=(20,20))
m = Basemap(llcrnrlat=24,urcrnrlat=50,\
            llcrnrlon=-125.5,urcrnrlon=-66.5,lat_ts=50,resolution='i',area_thresh=10000)
m.drawcoastlines()
m.drawcountries()
m.drawstates()
m.drawcounties()
x,y=m(nearestStationLon,nearestStationLat)
plotHandle = m.scatter(x,y,c='r',zorder=10)    
savefig('figures/nearestWeatherStation')
show()

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################

import pickle
from pylab import *
from sys import exit
import numpy as np
from datetime import datetime
from datetime import timedelta

latall=pickle.load(open('pickle_files/cmiplatall.p','rb'))
lonall=pickle.load(open('pickle_files/cmiplonall.p','rb'))

modelData85_tmax=pickle.load(open('pickle_files/modelData_all_tmax.p','rb'))
modelData85_tmin=pickle.load(open('pickle_files/modelData_all_tmin.p','rb'))

modelData45_tmax=pickle.load(open('pickle_files/modelData_all_tmax_rcp45.p','rb'))
modelData45_tmin=pickle.load(open('pickle_files/modelData_all_tmin_rcp45.p','rb'))

countyName=pickle.load(open('pickle_files/countyName.p','rb'))
startGrowingMon=pickle.load(open('pickle_files/startGrowingMon.p','rb'))
startGrowingDay=pickle.load(open('pickle_files/startGrowingDay.p','rb'))
endGrowingMon=pickle.load(open('pickle_files/endGrowingMon.p','rb'))
endGrowingDay=pickle.load(open('pickle_files/endGrowingDay.p','rb'))
cIndex=pickle.load(open('pickle_files/cIndex.p','rb'))
cIndexState=pickle.load(open('pickle_files/cIndexState.p','rb'))
T90all=pickle.load(open('pickle_files/T90all.p','rb'))
T10all=pickle.load(open('pickle_files/T10all.p','rb'))
slope=pickle.load(open('pickle_files/slope.p','rb'))
bIntercept=pickle.load(open('pickle_files/bIntercept.p','rb'))
presentGrowingCounties=pickle.load(open('pickle_files/presentGrowingCounties.p','rb'))
timeday=np.arange(42368,73049)
cropTitle=('corn','soybeans','rice')
runCrops=zeros(shape=(3))
icity=-1
nCities=3143 # number of cities  
nCrop=3
nyears=100
nScen=2
nPredictor=3

c=0

futureYield=-9999*ones(shape=(nCities,nPredictor,nyears,nScen,nCrop))

HeatWaves=-9999*ones(shape=(nScen,3143,100))
SummerAvg=-9999*ones(shape=(nScen,3143,100))
KDDays=-9999*ones(shape=(nScen,3143,100))

for i in range(nCities):
    icity+=1
    
    if icity!=604:
        continue
        
    if sum(presentGrowingCounties[icity])==0:
        continue
        
    doCorn=1
    doSoy=1
    doRice=1
    
    if presentGrowingCounties[icity,0]==0:
        doCorn=0
    if presentGrowingCounties[icity,1]==0:
        doSoy=0
    if presentGrowingCounties[icity,2]==0:
        doRice=0
        
    runCrops[0]=doCorn
    runCrops[1]=doSoy
    runCrops[2]=doRice
        
    lat=latall[icity]
    lon=lonall[icity]
    
    if int(lat)==-9999:
        c+=1
        continue
    
    city=countyName[icity].title()
    iBeg=16
    iEnd=100

    print icity,'of',nCities; sys.stdout.flush()

    dailyData=-9999*ones(shape=(nScen,2,nyears,12,31)) #v,y,m,d,scen
    
    for iday in range(len(timeday)):
        date_format = "%m/%d/%Y"
        d = timedelta(days=int(timeday[iday]))
        a = datetime.strptime('1/1/1900', date_format)
        date=a+d
        year=date.year
        month=date.month
        day=date.day
        
        y=year-2000
        y2=year-2016
        m=month-1
        d=day-1
           
        dailyData[0,0,y,m,d]=modelData85_tmax[d,m,y2,icity]
        dailyData[0,1,y,m,d]=modelData85_tmin[d,m,y2,icity]
        
        dailyData[1,0,y,m,d]=modelData45_tmax[d,m,y2,icity]
        dailyData[1,1,y,m,d]=modelData45_tmin[d,m,y2,icity]
    
    ###############################################
    # Make averages
    ###############################################
    mAvg=zeros(shape=(nScen,nyears,12))  #monthly average
    v=0
    
    ## compute yearly averages ##
    for scen in range(nScen):
        for y in range(iBeg,iEnd):
            for m in range(12):
                j=0
                for d in range(31):
                    if dailyData[scen,v,y,m,d]!=-9999:
                        mAvg[scen,y,m]+=dailyData[scen,v,y,m,d]
                        j+=1
                mAvg[scen,y,m]=mAvg[scen,y,m]/j
    
            
    yAvg=zeros(shape=(nScen,nyears))  #yearly average
    sAvg=zeros(shape=(nScen,nyears,4))  #seasonal average
    gAvg=zeros(shape=(nScen,nyears,nCrop)) #growing season average
    
    ## compute seasonal averages ##
    for scen in range(nScen):
        for y in range(iBeg,iEnd):
            yAvg[scen,y]=sum(mAvg[scen,y,:])/12
            
            sAvg[scen,y,0]=sum(mAvg[scen,y,2:5])/3 #spring MAM
            sAvg[scen,y,1]=sum(mAvg[scen,y,5:8])/3 #summer JJA
            sAvg[scen,y,2]=sum(mAvg[scen,y,8:11])/3 #fall SON
            sAvg[scen,y,3]=sum(mAvg[scen,y,0:2]+mAvg[scen,y,11])/3 #winter DJF
            
            SummerAvg[scen,icity,y]=sAvg[scen,y,1]
            
    ###############################################
    # Find Extremes
    ###############################################
    unsortedDall=zeros(shape=(nScen,4,366*315))
    # initialize variables for the extremes
    T90count=zeros(shape=(nScen,2,nyears))
    heatwavecount=zeros(shape=(nScen,2,nyears))
    seasonheatwavecount=zeros(shape=(nScen,2,nyears,nCrop))
    
    ## compute the extremes ##
    for scen in range(2):
        for v in range(1):
            i=3
            daysSinceHeatWave=0
            daysSinceColdWave=0
            for y in range(iBeg,iEnd):
                for m in range(12):
                    for d in range(31):
                        if dailyData[scen,v,y,m,d]>-100:
                            #unsortedDall[scen,v,i]=dailyData[scen,v,y,m,d]
                            unsortedDall[scen,v,i]=150
                            i+=1
                            daysSinceHeatWave+=1
                
                            if min(unsortedDall[scen,v,i-3:i])>=T90all[icity,v] and daysSinceHeatWave>2:
                                daysSinceHeatWave>2
                                heatwavecount[scen,v,y]+=1
                                daysSinceHeatWave=0
                                for cp in range(nCrop):
                                    if runCrops[cp]==0:
                                        continue
                                    if m>startGrowingMon[cIndexState[icity],cp] and m<endGrowingMon[cIndexState[icity],cp]:
                                        seasonheatwavecount[scen,v,y,cp]+=1
                                    
                                    elif m==startGrowingMon[cIndexState[icity],cp] and d>=startGrowingDay[cIndexState[icity],cp]:
                                        seasonheatwavecount[scen,v,y,cp]+=1
                                        
                                    elif m==endGrowingMon[cIndexState[icity],cp] and d<=endGrowingDay[cIndexState[icity],cp]:
                                        seasonheatwavecount[scen,v,y,cp]+=1
    
                HeatWaves[scen,icity,y]=heatwavecount[scen,0,y]
                                 
    ###############################################
    # Find KDD
    ###############################################  
    Tavg=zeros(shape=(nScen,nyears,12,31))
    KDD=zeros(shape=(nScen,nyears,nCrop))
    
    for scen in range(nScen):    
        for y in range(iBeg,iEnd):
            for m in range(12):
                for d in range(31):
                    if dailyData[scen,0,y,m,d]==-9999 or dailyData[scen,1,y,m,d]==-9999:
                        continue
                    dailytmin=dailyData[scen,1,y,m,d]
                    if dailytmin<50:
                        dailytmin=50
                    Tavg[scen,y,m,d]=(dailyData[scen,0,y,m,d]+dailytmin)/2
                    if Tavg[scen,y,m,d]<68:
                        continue
                    
                    for cp in range(nCrop):
                        if runCrops[cp]==0:
                                        continue
                        if m>startGrowingMon[cIndexState[icity],cp] and m<endGrowingMon[cIndexState[icity],cp]:
                            KDD[scen,y,cp]+=Tavg[scen,y,m,d]-68
                        
                        if m==startGrowingMon[cIndexState[icity],cp] and d>=startGrowingDay[cIndexState[icity],cp]:
                            KDD[scen,y,cp]+=Tavg[scen,y,m,d]-68  
                            
                        if m==endGrowingMon[cIndexState[icity],cp] and d<=endGrowingDay[cIndexState[icity],cp]:
                            KDD[scen,y,cp]+=Tavg[scen,y,m,d]-68 
    
            KDDays[scen,icity,y]=KDD[scen,y,0]
                           
    ###############################################
    # Predict Yields
    ###############################################
    ipredict=0    
    
    #print '\nKDD prediction:'
    for scen in range(nScen):
        #print 'scen=',scen
        for y in range(iBeg,iEnd):
            for cp in range(nCrop):
                if runCrops[cp]==0:
                    continue
                #print cropTitle[cp]   
                ## KDD prediction ##
                x=KDD[scen,:,cp]
                futureYield[icity,ipredict,y,scen,cp]=slope[2,icity,cp]*x[y]+bIntercept[2,icity,cp]
                futureYield[icity,ipredict,y,scen,cp]=round(futureYield[icity,ipredict,y,scen,cp],2)
                if -20<futureYield[icity,ipredict,y,scen,cp]<0:
                    futureYield[icity,ipredict,y,scen,cp]=0.01
                #print futureYield[icity,ipredict,y,scen,cp]
        #print '\n\n'
    ipredict+=1
    
    #print '\n\n\nSummer avg prediction:'
    for scen in range(nScen):
        #print 'scen=',scen
        for y in range(iBeg,iEnd):
            for cp in range(nCrop):
                if runCrops[cp]==0:
                    continue
                #print cropTitle[cp]
                    
                ## Summer avg prediction ##
                x=sAvg[scen,:,1]
                futureYield[icity,ipredict,y,scen,cp]=slope[3,icity,cp]*x[y]+bIntercept[3,icity,cp]
                futureYield[icity,ipredict,y,scen,cp]=round(futureYield[icity,ipredict,y,scen,cp],2)
                if -20<futureYield[icity,ipredict,y,scen,cp]<0:
                    futureYield[icity,ipredict,y,scen,cp]=0.01
                #print futureYield[icity,ipredict,y,scen,cp]
        #print '\n\n'
    ipredict+=1
    
    #print '\n\n\nHeat Wave prediction:'
    for scen in range(nScen):
        #print 'scen=',scen
        for y in range(iBeg,iEnd):
            for cp in range(nCrop):
                if runCrops[cp]==0:
                    continue
                #print cropTitle[cp]
                    
                ## Heat Wave prediction ##
                x=seasonheatwavecount[scen,0,:,cp]
                futureYield[icity,ipredict,y,scen,cp]=slope[4,icity,cp]*x[y]+bIntercept[4,icity,cp]
                futureYield[icity,ipredict,y,scen,cp]=round(futureYield[icity,ipredict,y,scen,cp],2)
                if -20<futureYield[icity,ipredict,y,scen,cp]<0:
                    futureYield[icity,ipredict,y,scen,cp]=0.01
                #print futureYield[icity,ipredict,y,scen,cp]
        #print '\n\n'

        
pickle.dump(futureYield,open('pickle_files/futureYield170105.p','wb'))  
pickle.dump(SummerAvg,open('pickle_files/SummerAvgFuture.p','wb'))
pickle.dump(HeatWaves,open('pickle_files/HeatWavesFuture.p','wb'))
pickle.dump(KDDays,open('pickle_files/KDDFuture.p','wb'))    







#imshow(tmax2[0,:,:],origin='lower',extent=[lon[0],lon[540],lat[0],lat[382]],vmin=255,vmax=300)
#colorbar()
#savefig('final_figures/model_data')
#show()

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################
 
################################################
# Makes the aggregated plots
# Makes the slope plots
################################################
from pylab import *
import csv
import math
from sys import exit
import pickle
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from mpl_toolkits.basemap import Basemap
import os
from scipy.stats import norm
import matplotlib as mpl
from matplotlib.patches import Polygon

###############################################
# Functions
###############################################
def reverse_colormap(cmap, name = 'my_cmap_r'):
    reverse = []
    k = []  
    for key in cmap._segmentdata:    
        k.append(key)
        channel = cmap._segmentdata[key]
        data = []

        for t in channel:                    
            data.append((1-t[0],t[2],t[1]))            
        reverse.append(sorted(data))    

    LinearL = dict(zip(k,reverse))
    my_cmap_r = mpl.colors.LinearSegmentedColormap(name, LinearL) 
    return my_cmap_r 
    
###############################################
# Pickle Files
###############################################    
cityFile='tmax_cities_over_45_years'
cIndexState=pickle.load(open('pickle_files/cIndexState.p','rb'))
Corr=pickle.load(open('pickle_files/Corr.p','rb'))
maxCorr=pickle.load(open('pickle_files/maxCorr.p','rb'))
iplotCorr=pickle.load(open('pickle_files/iplotCorr.p','rb'))
maxCorr2=pickle.load(open('pickle_files/maxCorr2.p','rb'))
iplotCorr2=pickle.load(open('pickle_files/iplotCorr2.p','rb'))
slope=pickle.load(open('pickle_files/slope.p','rb'))
R_squared=pickle.load(open('pickle_files/'+cityFile+'/R2.p','rb'))
lat=pickle.load(open('pickle_files/'+cityFile+'/lat.p','rb'))
lon=pickle.load(open('pickle_files/'+cityFile+'/lon.p','rb'))
station=pickle.load(open('pickle_files/'+cityFile+'/station.p','rb'))
goodCity=pickle.load(open('pickle_files/'+cityFile+'/goodCity.p','rb'))
cIndex=pickle.load(open('pickle_files/cIndex.p','rb'))
presentGrowingCounties=pickle.load(open('pickle_files/presentGrowingCounties.p','rb'))
highTotalYield=pickle.load(open('pickle_files/highTotalYield.p','rb'))
goodCity3=pickle.load(open('pickle_files/goodCity3.p','rb'))

###############################################
# Define Variables
###############################################
v=0

lat=lat[goodCity[v,:]]
lon=lon[goodCity[v,:]]
numGoodCities=sum(goodCity[v,:])
dataYield=ones(shape=(3143,2),dtype=bool)
cropTitle=('Corn','Soybeans','Rice')

#var=('Time','Normalized_Time','Growing_Degree_Days','Killing_Degree_Days','Frost_Nights',
#    'Season_Frost_Nights','Tmax_Avg','Tmin_Avg','Spring_Avg','Summer_Avg','Tropical_Nights',
#    'Tmax_Heat_Waves','Tmin_Heat_Waves','Tmax_Warm_Days','Tmin_Warm_Days')

var=('Time','Normalized_Time','Killing_Degree_Days','Summer_Avg','Heat_Waves')
titleVar=('Time','Normalized Time','Killing Degree Days','Summer Avg Temp','Heat Waves')
   
cp=0
cmapArray=plt.cm.hot(arange(256))
cmin=-0.7
cmax=0
y1=0
y2=255

MinMaxArray=ones(shape=(2,3))

b=0
significant=0
highlySignificant=0
countyCounter=0
###############################################
# Plot
###############################################
MinMaxArray=zeros(shape=(2,1))
goodCity2=zeros(shape=(3143,3),dtype=bool)

for iplot in range(4,5):
    print iplot
    
    figure(5,figsize=(9,7))
    show()
    
    for icity in range(3143):
        if Corr[iplot,icity,cp]!=0:
            goodCity2[icity,cp]=True
    
    j=-1
    
    subPlot1 = plt.axes([.5,.25,.475,.5])
    MinMaxArray[0,0]=-0.7
    MinMaxArray[1,0]=0
    
    plt.imshow(MinMaxArray,cmap='hot')
    plt.colorbar()
    
    # create the map
    subPlot1 = plt.axes([0.1,0,0.8,1])
    m = Basemap(llcrnrlon=-107,llcrnrlat=27,urcrnrlon=-71,urcrnrlat=49,
    lat_ts=50,resolution='i',area_thresh=10000)
            
    # load the shapefile, use the name 'states'
    m.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
    
    ax = plt.gca() # get current axes instance
    
    for shape_dict in m.states_info:
        j+=1
        seg = m.states[j]
        
        s=int(shape_dict['STATEFP'])
        c=int(shape_dict['COUNTYFP'])
        
        if s==72 or cIndex[s,c]==-9999:
            continue
            
        if goodCity2[cIndex[s,c],cp]==False:
            continue
            
        x=Corr[iplot,cIndex[s,c],cp]
        if math.isnan(x)==True:
            b+=1
            goodCity2[cIndex[s,c],cp]=False
            continue
        countyCounter+=1
        if x<=-.49:
            significant+=1
        if x<=-.59:
            highlySignificant+=1
        y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
        icmap=min(255,int(round(y,1)))
        
         
        poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
        ax.add_patch(poly)
    #print titleVar[iplot],round(mean(Corr[iplot,goodCity2[:,cp],cp]),2),significant, round((float(significant)/75.)*100,2)
    #print titleVar[iplot],round(mean(Corr[iplot,goodCity2[:,cp],cp]),2),highlySignificant, round((float(highlySignificant)/75.)*100,2)
    title('Correlation of '+cropTitle[cp]+' and '+titleVar[iplot]+', Avg Corr= '+str(round(mean(Corr[iplot,goodCity2[:,cp],cp]),2)))
    savefig('final_figures/correlation/'+cropTitle[cp]+'/'+str(var[iplot])+'_'+cropTitle[cp]+'_yield_corr',dpi=500)
    show() 
    clf()

    
           
###############################################
# Plot Slopes
###############################################
cp=2

cmapArray=plt.cm.gist_rainbow(arange(256))
y1=0
y2=255

MinMaxArray=ones(shape=(2,3))

b=0
significant=0
highlySignificant=0
countyCounter=0

MinMaxArray=zeros(shape=(2,1))
goodCity2=zeros(shape=(3143,3),dtype=bool)

for iplot in range(3,5):
    print iplot
    
    figure(5,figsize=(9,7))
    show()
    
    for icity in range(3143):
        if slope[iplot,icity,cp]!=0 or slope[iplot,icity,cp]!=-9999:
            goodCity2[icity,cp]=True
    
    cmax=np.amax(slope[iplot,highTotalYield[:,cp],cp])
    #cmax=0
    cmin=np.amin(slope[iplot,highTotalYield[:,cp],cp])
    if cp==0:
        if iplot==2:
            cmin=-0.1*10
            cmax=0.02*10
        if iplot==3:
            cmin=-8
            cmax=1.25
        if iplot==4:
            cmin=-4
            cmax=.7
            
    if cp==1:
        if iplot==2:
            cmin=-0.02*10
            cmax=0.006*10
        if iplot==3:
            cmin=-1.3
            cmax=.3
        if iplot==4:
            cmin=-.8
            cmax=.2
            
    if cp==2:
        if iplot==3:
            cmax=20
            cmin=-85
        if iplot==4:
            cmax=15
            cmin=-50
            
    if cmin==-9999:
        print 'cmin=-9999'
    
    j=-1
    
    subPlot1 = plt.axes([.5,.25,.475,.5])
    MinMaxArray[0,0]=cmax
    MinMaxArray[1,0]=cmin
    
    plt.imshow(MinMaxArray,cmap='gist_rainbow')
    plt.colorbar()
    
    # create the map
    subPlot1 = plt.axes([0.1,0,0.8,1])
    m = Basemap(llcrnrlon=-107,llcrnrlat=27,urcrnrlon=-71,urcrnrlat=49,
    lat_ts=50,resolution='i',area_thresh=10000)
            
    # load the shapefile, use the name 'states'
    m.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
    
    ax = plt.gca() # get current axes instance
    
    for shape_dict in m.states_info:
        j+=1
        seg = m.states[j]
        
        s=int(shape_dict['STATEFP'])
        c=int(shape_dict['COUNTYFP'])
        
        if s==72 or cIndex[s,c]==-9999:
            continue
            
        if highTotalYield[cIndex[s,c],cp]==False:
            continue
            
        x=slope[iplot,cIndex[s,c],cp]
        if math.isnan(x)==True:
            b+=1
            goodCity2[cIndex[s,c],cp]=False
            continue

        y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
        icmap1=min(255,int(round(y,1)))
        icmap=max(0,icmap1)
        
         
        poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
        ax.add_patch(poly)
    title('Slope: '+cropTitle[cp]+' Yield/'+titleVar[iplot]+', Avg = '+str(round(mean(slope[iplot,highTotalYield[:,cp],cp]),2)))
    savefig('final_figures/slope/'+cropTitle[cp]+'/'+str(var[iplot])+'_'+cropTitle[cp]+'_yield_slope',dpi=500)
    show() 
    pause(1)
    clf()

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################
   
################################################
# Make movie of Yield summer avg over time
################################################
from pylab import *
import csv
from math import sqrt
from sys import exit
import pickle
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from mpl_toolkits.basemap import Basemap
import os
from scipy.stats import norm
import matplotlib as mpl
from matplotlib.patches import Polygon

###############################################
# Functions
###############################################
def Avg(x):   
    '''function to average'''
    xAvg=0.
    for k in range(len(x)):
        xAvg=xAvg+x[k]
    xAvg=xAvg/(k+1)
    return xAvg 
    
    
###############################################
# Variables
###############################################
futureYield=pickle.load(open('pickle_files/futureYield.p','rb'))
# futureYield dimensions = (nCities,nPredictor,nyears,nScen,nCrop)
# nPredictor  0=Summer avg, 1=Heat Waves, 2=KDD
cropYield=pickle.load(open('pickle_files/cropYield.p','rb'))
# cropYield dimensions = (nCities,year,cp)
cIndex=pickle.load(open('pickle_files/cIndex.p','rb'))
cIndexState=pickle.load(open('pickle_files/cIndexState.p','rb'))
BBcounty=pickle.load(open('pickle_files/BBcounty.p','rb'))
SummerAvg=pickle.load(open('pickle_files/SummerAvg.p','rb'))
KDD=pickle.load(open('pickle_files/KDD.p','rb'))
HeatWaves=pickle.load(open('pickle_files/HeatWaves.p','rb'))
presentGrowingCounties=pickle.load(open('pickle_files/presentGrowingCounties.p','rb'))
KDDFuture=pickle.load(open('pickle_files/KDDFuture.p','rb'))
HeatWavesFuture=pickle.load(open('pickle_files/HeatWavesFuture.p','rb'))
SummerAvgFuture=pickle.load(open('pickle_files/SummerAvgFuture.p','rb'))
highTotalYield=pickle.load(open('pickle_files/highTotalYield.p','rb'))

goodfuture=zeros(shape=(3143,3),dtype=bool)
goodpast=zeros(shape=(3143,3),dtype=bool)
goodpastyears=zeros(shape=(3143,116,3),dtype=bool)

avgFutureYield=zeros(shape=(100,3))
avgCropYield=zeros(shape=(116,3)) 
nCrop=3  
cropTitle=('Corn','Soybeans','Rice')

for icity in range(3143):
    for y in range(16,100):
        for cp in range(nCrop):
            if futureYield[icity,0,y,0,cp]>0 and presentGrowingCounties[icity,cp]==1 and isnan(futureYield[icity,0,y,0,cp])==False:
                goodfuture[icity,cp]=True
    for y in range(70,116):
        for cp in range(nCrop):
            if cropYield[icity,y,cp]>0 and highTotalYield[icity,cp]==True:
                    goodpast[icity,cp]=True
                    goodpastyears[icity,y,cp]=True
       
###############################################
# Plot Yield 1970-2015
###############################################
for cp in range(3):
    cmapArray=plt.cm.jet(arange(256))
    cmin=0
    y1=0
    y2=256
    year=115
    cmax=.8*(np.amax(cropYield[:,105:116,cp]))
    
    for year in range(70,116):
        print year+1900
        
        figure(1,figsize=(9,7))
        show()
        
        map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
                projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
        
        # load the shapefile, use the name 'states'
        map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
        
        ax = plt.gca() # get current axes instance
        j=-1
        
        for shape_dict in map.states_info:
            j+=1
            seg = map.states[j]
            
            s=int(shape_dict['STATEFP'])
            c=int(shape_dict['COUNTYFP'])
            
            if s==72 or cIndex[s,c]==-9999:
                continue
        
            Yield=cropYield[cIndex[s,c],year,cp]
            
            if Yield<=0:
                continue
                
            x=Yield
            y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
            icmap=min(255,int(round(y,1)))
            
            poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
            ax.add_patch(poly)
            
        title('Bushels/Acre of '+cropTitle[cp]+' Grown')
        text(3500000,2500000,year+1900,fontsize=25)
        savefig('final_figures/crop_yield_every_year/'+cropTitle[cp]+'/'+cropTitle[cp]+'_yield_'+str(year+1900))
        plt.show()      
        clf()
    
###############################################
# Plot Yield 2016-2100
###############################################
for cp in range(nCrop):
    cmapArray=plt.cm.jet(arange(256))
    cmin=0
    y1=0
    y2=256
    cmax=.8*(np.amax(cropYield[:,105:116,cp]))
    
    for year in range(16,100):
        print year+2000
        
        figure(1,figsize=(9,7))
        show()
        
        map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
                projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
        
        # load the shapefile, use the name 'states'
        map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
        
        ax = plt.gca() # get current axes instance
        j=-1
        
        for shape_dict in map.states_info:
            j+=1
            seg = map.states[j]
            
            s=int(shape_dict['STATEFP'])
            c=int(shape_dict['COUNTYFP'])
            
            if s==72 or cIndex[s,c]==-9999:
                continue
        
            Yield=mean(futureYield[cIndex[s,c],0:2,year,0,cp])
            
            if Yield<=0 or isnan(Yield)==True:
                continue
                
            x=Yield
            y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
            icmap=min(255,int(round(y,1)))
            
            poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
            ax.add_patch(poly)
            
        title('Bushels/Acre of '+cropTitle[cp]+' Grown in '+str(year+2000))
        text(3500000,2500000,year+2000,fontsize=25)
        savefig('final_figures/crop_yield_every_year/'+cropTitle[cp]+'/'+cropTitle[cp]+'_yield_'+str(year+2000))
        plt.show() 
        clf()

###############################################
# Plot Summer Avg Temp
###############################################
cmapArray=plt.cm.jet(arange(256))
cmin=74
y1=0
y2=256
cmax=98

for year in range(70,115):
    b=0
    print year+1900
    
    figure(1,figsize=(9,7))
    show()
    
    map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
    
    # load the shapefile, use the name 'states'
    map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
    
    ax = plt.gca() # get current axes instance
    j=-1
    
    for shape_dict in map.states_info:
        j+=1
        seg = map.states[j]
        
        s=int(shape_dict['STATEFP'])
        c=int(shape_dict['COUNTYFP'])
        
        if s==72 or cIndex[s,c]==-9999:
            b+=1
            continue
    
        Temp=SummerAvg[cIndex[s,c],year]
        
        if Temp<=0:
            b+=1
            continue
            
        x=Temp
        y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
        icmap=min(255,int(round(y,1)))
        
        poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
        ax.add_patch(poly)
        
    title('Summer Average Temperature')
    text(3500000,2500000,year+1900,fontsize=25)
    savefig('final_figures/summer_avg_every_year/summer_'+str(year+1900))
    plt.show() 
    clf()

###############################################
# Plot Summer Avg Temp 2016-2100
###############################################
cmapArray=plt.cm.jet(arange(256))
cmin=74
y1=0
y2=256
cmax=98

for year in range(16,100):
    b=0
    print year+2000
    
    figure(1,figsize=(9,7))
    show()
    
    map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
    
    # load the shapefile, use the name 'states'
    map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
    
    ax = plt.gca() # get current axes instance
    j=-1
    
    for shape_dict in map.states_info:
        j+=1
        seg = map.states[j]
        
        s=int(shape_dict['STATEFP'])
        c=int(shape_dict['COUNTYFP'])
        
        if s==72 or cIndex[s,c]==-9999:
            b+=1
            continue
    
        Temp=SummerAvgFuture[cIndex[s,c],year]
        
        if Temp<=0 or isnan(Temp)==True:
            b+=1
            continue
            
        x=Temp
        y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
        icmap=min(255,int(round(y,1)))
        
        poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
        ax.add_patch(poly)
        
    title('Summer Average Temperature in '+str(year+2000))
    text(3500000,2500000,year+2000,fontsize=25)
    savefig('final_figures/summer_avg_every_year/summer_'+str(year+2000))
    plt.show() 
    clf()

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################

###############################################
# The Main Code
# Reads in data
# gets rid of bad cities
# computes the means and extremes
# plots and records slope, corr, lat, lon, station name
###############################################
from pylab import *
import csv
from math import sqrt
from sys import exit
import numpy as np
import pickle
import os
from scipy import stats
import matplotlib.lines as lines

# cd Documents/Science_Fair_2017_Crop_Yields/

# choose what you want the code to run for #
useUniformStartDate=True
dataYears=45
goodCityRatio=.9
makePlots=True
nCrop=3

###############################################
# Functions
###############################################
def Avg(x):   
    '''function to average'''
    xAvg=0.
    for k in range(len(x)):
        xAvg=xAvg+x[k]
    xAvg=xAvg/(k+1)
    return xAvg 
    
def stdDev(x):   
    '''function to compute standard deviation'''
    xAvg=Avg(x)
    xOut=0.
    for k in range(len(x)):
        xOut=xOut+(x[k]-xAvg)**2
    xOut=xOut/(k+1)
    xOut=sqrt(xOut)
    return xOut
    
def Variance(x):   
    '''function to compute the variance (std dev squared)'''
    xAvg=Avg(x)
    xOut=0.
    for k in range(len(x)):
        xOut=xOut+(x[k]-xAvg)**2
    xOut=xOut/(k+1)
    return xOut
    
def SumOfSquares(x):
    '''function to compute the sum of squares'''
    xOut=0.
    for k in range(len(x)):
        xOut=xOut+x[k]**2
    return xOut

def corr(x,y):   
    ''' function to find the correlation of two arrays'''
    xAvg=Avg(x)
    Avgy=Avg(y)
    rxy=0.
    n=min(len(x),len(y))
    for k in range(n):
        rxy=rxy+(x[k]-xAvg)*(y[k]-Avgy)
    rxy=rxy/(k+1)
    stdDevx=stdDev(x)
    stdDevy=stdDev(y)
    rxy=rxy/(stdDevx*stdDevy)
    return rxy
    
###############################################
# Read in Weather Data
###############################################
cityFile='tmax_cities_over_'+str(dataYears)+'_years' # what file to read in #

nCities=3143 # number of cities   
#f=open('written_files/'+cityFile+'.txt','r')
f=open('written_files/tmax_nearest_stations.txt')

# read in the pickle files
startGrowingMon=pickle.load(open('pickle_files/startGrowingMon.p','rb'))
startGrowingDay=pickle.load(open('pickle_files/startGrowingDay.p','rb'))
endGrowingMon=pickle.load(open('pickle_files/endGrowingMon.p','rb'))
endGrowingDay=pickle.load(open('pickle_files/endGrowingDay.p','rb'))
cIndex=pickle.load(open('pickle_files/cIndex.p','rb'))
cIndexState=pickle.load(open('pickle_files/cIndexState.p','rb'))
cropYield=pickle.load(open('pickle_files/cropYield.p','rb'))
countyName=pickle.load(open('pickle_files/countyName.p','rb'))
presentGrowingCounties=pickle.load(open('pickle_files/presentGrowingCounties.p','rb'))

# initialize variables
icity=-1 # a counter of the number of cities
slope=-9999.*ones(shape=(5,nCities,nCrop))
bIntercept=-9999.*ones(shape=(5,nCities,nCrop))
lat=-9999.*ones(shape=(nCities))
lon=-9999.*ones(shape=(nCities))
lat2=-9999.*ones(shape=(nCities))
lon2=-9999.*ones(shape=(nCities))
stationList=[]
cityList=[]
goodCity=ones(shape=(4,nCities),dtype=bool)
Corr=zeros(shape=(5,nCities,nCrop))
R2=-9999.*ones(shape=(5,nCities,nCrop))
maxCorr=zeros(shape=(nCities,nCrop))
iplotCorr=zeros(shape=(nCities,nCrop))
maxCorr2=zeros(shape=(nCities,nCrop))
iplotCorr2=zeros(shape=(nCities,nCrop))
T90all=zeros(shape=(nCities,2))
T10all=zeros(shape=(nCities,2))
runCrops=zeros(shape=(3))
dataAt2015=-9999*ones(shape=(nCities,nCrop))

HeatWaves=-9999*ones(shape=(3143,116))
SummerAvg=-9999*ones(shape=(3143,116))
KDDays=-9999*ones(shape=(3143,116))
mostcity=[600,604,632,644,647,651,843]

for cityline in f:  # read in the file with the closest counties
    icity+=1
    
    if sum(presentGrowingCounties[icity])==0:
        continue
    
    
    if icity!=604:
        continue
        
    doCorn=1
    doSoy=1
    doRice=1
        
    # initialize more variables
    cityline=cityline.translate(None, ' ')
    tmp=cityline.split(',')
    station=tmp[0]
    stationList.append(station)
    lat[icity]=tmp[2]
    lon[icity]=tmp[3]
    iBeg=float(tmp[5])-1900
    iBeg=int(iBeg)
    if useUniformStartDate:
        iBeg=(2015-dataYears)-1900
    iEnd1=float(tmp[6])-1900
    iEnd1=int(iEnd1)
    city=tmp[7].strip()
    city=countyName[icity].title()
     
    station2=tmp[8]
    lat2[icity]=tmp[10]
    lon2[icity]=tmp[11]
    iEnd2=float(tmp[14])-1900
    iEnd2=int(iEnd2)

    iEnd=max(iEnd1,iEnd2)
    
    print '\n'
    print icity,'of',nCities
    print city
   
    # initialize more variables
    dailyData1=-9999.*ones(shape=(4,116,12,31))
    badData=zeros(shape=(4,116))  
    badYears=ones(shape=(4,116))
    badYearsforAvg=ones(shape=(4,116,nCrop))
    badYearYield=ones(shape=(4,116,nCrop))
    badYearYieldBoolean=zeros(shape=(4,116,nCrop),dtype=bool)
    badYear=ones(shape=(4,116))
        
    dailyData2=-9999.*ones(shape=(4,116,12,31))
    badData2=zeros(shape=(4,116))  
        
    f = open('data/noaa_daily_data/ghcnd_all/'+station+'.dly', 'r')
    f2=open('data/noaa_daily_data/ghcnd_all/'+station2+'.dly', 'r')
    
    for line in f: #read in the daily data for each closest station
        # initialize tmp variables
        tmp=line[0:270]
        station=tmp[0:11]
        year=float(tmp[11:15])
        month=float(tmp[15:17])
        var=tmp[17:21]
        if var=='TMAX':
            v=0
        elif var=='TMIN':
            v=1
        else:
            v=-1
        
        if var!=-1:   
            m=month-1
            y=year-1900
            d=21
            if y<0:
                y=0
            for d in range(31):
                tmp2=tmp[21+8*d:21+8*d+5]
                if tmp2=='   0T' or tmp2=='   0P':
                    tmp2='    0'
                if tmp2!='-9999':
                    if v==0 or v==1:
                        tmp3=(9./5.*float(tmp2)/10.+32.)
                    else:
                        tmp3=(float(tmp2)/10/10/2.54)
                else:
                    badData[v,y]+=1
                    tmp3=-9999
                dailyData1[v,y,m,d]=tmp3 # put the data into one 4D array
    
    
    for line in f2: #read in the daily data for each second closest station
        # initialize tmp variables
        tmp=line[0:270]
        station=tmp[0:11]
        year=float(tmp[11:15])
        month=float(tmp[15:17])
        var=tmp[17:21]
        if var=='TMAX':
            v=0
        elif var=='TMIN':
            v=1
        elif var=='PRCP':
            v=2
        elif var=='SNOW':
            v=3
        else:
            v=-1
        
        if var!=-1:   
            m=month-1
            y=year-1900
            d=21
            if y<0:
                y=0
            for d in range(31):
                tmp2=tmp[21+8*d:21+8*d+5]
                if tmp2=='   0T' or tmp2=='   0P':
                    tmp2='    0'
                if tmp2!='-9999':
                    if v==0 or v==1:
                        tmp3=(9./5.*float(tmp2)/10.+32.)
                    else:
                        tmp3=(float(tmp2)/10/10/2.54)
                else:
                    badData2[v,y]+=1
                    tmp3=-9999
                dailyData2[v,y,m,d]=tmp3 # put the data into one 4D array
    
    #################################################
    # Average the closest and second closest station
    #################################################
    dailyData=-9999*ones(shape=(4,116,12,31))
    mAvg=zeros(shape=(4,116,12))  #monthly average
    goodMonths=zeros(shape=(4,116,12))
    goodydcorn=0
    goodydsoy=0
    goodydrice=0
    
    for v in range(2):
        for y in range(iBeg,iEnd):
            g=0 # number of good months
            for m in range(12):
                j=0    # j is number of good days in month
                b=0    # b is the number of bad days in a month
                for d in range(31):
                    if dailyData1[v,y,m,d]>-100 and dailyData2[v,y,m,d]>-100: 
                        j+=1 
                        dailyData[v,y,m,d]=(dailyData1[v,y,m,d]+dailyData2[v,y,m,d])/2
                        
                    if dailyData1[v,y,m,d]>-100 and dailyData2[v,y,m,d]<-100:
                        j+=1
                        dailyData[v,y,m,d]=dailyData1[v,y,m,d]
                        
                    if dailyData1[v,y,m,d]<-100 and dailyData2[v,y,m,d]>-100:
                        j+=1
                        dailyData[v,y,m,d]=dailyData2[v,y,m,d]
                        
                    if dailyData1[v,y,m,d]<-100 and dailyData2[v,y,m,d]<-100 and d!=30:
                        b+=1
                        dailyData[v,y,m,d]=dailyData[v,y,m,d-1]
                if j>25:
                    g+=1
                        
                        
                    
#                if j>=20:
                    goodMonths[v,y,m]=1 # one means the month is good
                if j<20:   # make sure the month has atleast 20 good day
                    goodMonths[v,y,m]=0
                    mAvg[v,y,m]=-9999
            if goodMonths[v,y,3]==1 and goodMonths[v,y,4]==1 and goodMonths[v,y,5]==1 and goodMonths[v,y,6]==1 and goodMonths[v,y,7]==1 and goodMonths[v,y,8]==1 and goodMonths[v,y,9]==1 and goodMonths[v,y,10]==1:
                badYear[v,y]=0 # zero means the year is good

            if sum(goodMonths[v,y])==12:
                badYears[v,y]=0 # zero means the year is good
                
        if sum(badYear[v,iBeg:iEnd])/float(dataYears) > goodCityRatio: # if the city has too much bad data
            goodCity[v,icity]=False # True means the city is good
        
        for cp in range(nCrop):        
            for y in range(iBeg,iEnd):
                if cropYield[icity,y,cp]>0 and badYear[v,y]==0:
                    badYearYield[v,y,cp]=0 # zero means the year is good
                    badYearYieldBoolean[v,y,cp]=True
                if cropYield[icity,y,cp]!=-9999 and badYears[v,y]==0:
                    badYearsforAvg[v,y,cp]=0
                if cropYield[icity,y,cp]>0 and cp==0:
                    goodydcorn+=1
                if cropYield[icity,y,cp]>0 and cp==1:
                    goodydsoy+=1
                if cropYield[icity,y,cp]>0 and cp==2:
                    goodydrice+=1
                
           
    if goodCity[0,icity]==False:
        continue
    if goodydcorn<20:
        #print 'WARNING: NO CORN YIELD DATA'
        doCorn=0
    if goodydsoy<20:
        #print 'WARNING: NO SOY YIELD DATA'
        doSoy=0
    if goodydrice<20:
        #print 'WARNING: NO RICE YIELD DATA'
        doRice=0
    
    runCrops[0]=doCorn
    runCrops[1]=doSoy
    runCrops[2]=doRice
    ###############################################
    # Make averages
    ###############################################
    mAvg=zeros(shape=(4,116,12))  #monthly average
    year=range(1900,2016)
    maxP=-9999*ones(shape=(116))  # the maximum precip in one day for that year
    
    ## compute yearly averages ##
    for v in range(2):
        for y in range(iBeg,iEnd):
            for m in range(12):
                j=0             # j is number of good days in month
                for d in range(31):
                    if dailyData[v,y,m,d]>-100:
                        mAvg[v,y,m]+=dailyData[v,y,m,d]
                        j+=1   
                        if v==2:
                            maxP[y]=max(maxP[y],dailyData[v,y,m,d])
                if v==0 or v==1:
                    mAvg[v,y,m]=mAvg[v,y,m]/j
                if j==0:   # make sure the month has atleast one good day
                    mAvg[v,y,m]=-9999
    
            
    yAvg=zeros(shape=(4,116))  #yearly average
    sAvg=zeros(shape=(4,116,4))  #seasonal average
        
    ## compute seasonal averages ##
    for v in range(2):
        for y in range(iBeg,iEnd):
            if badYears[v,y]==0:
                if v==0 or v==1:
                    yAvg[v,y]+=sum(mAvg[v,y,:])/12
                else:
                    yAvg[v,y]+=sum(mAvg[v,y,:]) 
            
            if v==0 or v==1:
                sAvg[v,y,0]+=sum(mAvg[v,y,2:5])/3 #spring MAM
                sAvg[v,y,1]+=sum(mAvg[v,y,5:8])/3 #summer JJA
                sAvg[v,y,2]+=sum(mAvg[v,y,8:11])/3 #fall SON
                sAvg[v,y,3]+=sum(mAvg[v,y,0:2]+mAvg[v,y,11])/3 #winter DJF
            if sAvg[0,y,1]>0:
                SummerAvg[icity,y]=sAvg[0,y,1]
      
    ###############################################
    # Find percentiles
    ###############################################
    unsortedD=zeros(shape=(4,10970))
    sortedD=zeros(shape=(4,10970))
    
    
    for v in range(2):
        i=0
        for y in range(61,90):
            for m in range(12):
                for d in range(31):     
                    if dailyData[v,y,m,d]>-100:
                        unsortedD[v,i]=dailyData[v,y,m,d]
                        i+=1

        
        sortedD[v,0:i]=sort(unsortedD[v,0:i])
        i90=int(i*9/10.0)
        i10=int(i/10.)
        
    T90=array([sortedD[0,i90],sortedD[1,i90]]) # 90th percentile
    T10=array([sortedD[0,i10],sortedD[1,i10]]) # 10th percentile
    
    T90all[icity]=T90
    T10all[icity]=T10
    
    ###############################################
    # Find Extremes
    ###############################################
    year=range(1900,2016)
    nyears=max(year)-min(year)+1
    unsortedDall=zeros(shape=(4,366*315))
    # initialize variables for the extremes
    T90count=zeros(shape=(2,nyears))
    seasonT90count=zeros(shape=(2,nyears,nCrop))
    T10count=zeros(shape=(2,nyears))
    maxT=-9999*ones(shape=(2,nyears))
    minT=9999*ones(shape=(2,nyears))
    frostDays=zeros(shape=(nyears))
    seasonfrostDays=zeros(shape=(nyears,nCrop))
    tropicNights=zeros(shape=(nyears))
    seasontropicNights=zeros(shape=(nyears,nCrop))
    precip95count=zeros(shape=(nyears))
    heatwavecount=zeros(shape=(2,nyears))
    seasonheatwavecount=zeros(shape=(2,nyears,nCrop))
    coldwavecount=zeros(shape=(2,nyears))
    yr=arange(int(min(year)),int(max(year)+1))
    
    ## compute the extremes ##
    for v in range(2):
        i=3
        daysSinceHeatWave=0
        daysSinceColdWave=0
        for y in range(iBeg,iEnd):
            for m in range(12):
                for d in range(31):
                    if dailyData[v,y,m,d]>-100:
                        unsortedDall[v,i]=dailyData[v,y,m,d]
                        i+=1
                        daysSinceHeatWave+=1
                        daysSinceColdWave+=1    

                        if min(unsortedDall[v,i-3:i])>=T90[v] and daysSinceHeatWave>2:
                            daysSinceHeatWave>2
                            heatwavecount[v,y]+=1
                            daysSinceHeatWave=0
                            for cp in range(nCrop):
                                if runCrops[cp]==0:
                                    continue
                                if m>startGrowingMon[cIndexState[icity],cp] and m<endGrowingMon[cIndexState[icity],cp]:
                                    seasonheatwavecount[v,y,cp]+=1
                                
                                elif m==startGrowingMon[cIndexState[icity],cp] and d>=startGrowingDay[cIndexState[icity],cp]:
                                    seasonheatwavecount[v,y,cp]+=1
                                    
                                elif m==endGrowingMon[cIndexState[icity],cp] and d<=endGrowingDay[cIndexState[icity],cp]:
                                    seasonheatwavecount[v,y,cp]+=1
        
            if heatwavecount[0,y]>=0:
                HeatWaves[icity,y]=heatwavecount[0,y]
    ###############################################
    # Find GDD and KDD
    ###############################################  
    Tavg=zeros(shape=(116,12,31))
    GDD=zeros(shape=(116,nCrop))
    KDD=zeros(shape=(116,nCrop))
        
    for y in range(iBeg,iEnd):
        for m in range(12):
            for d in range(31):
                if dailyData[0,y,m,d]==-9999 or dailyData[1,y,m,d]==-9999:
                    continue
                dailytmin=dailyData[1,y,m,d]
                if dailytmin<50:
                    dailytmin=50
                Tavg[y,m,d]=(dailyData[0,y,m,d]+dailytmin)/2
                if Tavg[y,m,d]<68:
                    continue
                
                for cp in range(nCrop):
                    if runCrops[cp]==0:
                        continue
                    if m>startGrowingMon[cIndexState[icity],cp] and m<endGrowingMon[cIndexState[icity],cp]:
                        KDD[y,cp]+=Tavg[y,m,d]-68
                    
                    if m==startGrowingMon[cIndexState[icity],cp] and d>=startGrowingDay[cIndexState[icity],cp]:
                        KDD[y,cp]+=Tavg[y,m,d]-68  
                        
                    if m==endGrowingMon[cIndexState[icity],cp] and d<=endGrowingDay[cIndexState[icity],cp]:
                        KDD[y,cp]+=Tavg[y,m,d]-68                                    
                        
        if KDD[y,0]>=0:
            KDDays[icity,y]=KDD[y,0]                   
    ###############################################
    # Make Corr Plots
    ###############################################
    time=True           # iplot=0
    Normalizedtime=True # iplot=1  
    UseKDD=True         # iplot=2
    Summer=True         # iplot=3
    UseHeatWaves=True   # iplot=4
    
    DN=('Days','Nights')
    V=('Max Temp.','Min Temp.','Precip.','Snowfall')
    T=('Tmax','Tmin','precip','snow')
    HL=('Highs','Lows')
    capsDN=('DAYS','NIGHTS')
    capsDNsingle=('DAY','NIGHT')
    ncrop=('Corn','Soybean')
    iplot=0
    x=yr[iBeg:iEnd]
    if not os.path.exists(r'figures/mostcity/'+str(cIndexState[icity])+'/'+city):
        os.makedirs(r'figures/mostcity/'+str(cIndexState[icity])+'/'+city)
    
    ##                        ##
    ## Make Temperature Plots ##
    ##                        ##
    
    if time: # iplot=0
        for cp in range(nCrop):
            if runCrops[cp]==0:
                continue
            x=np.ma.compressed(np.ma.masked_array(yr[iBeg:iEnd],badYearYield[0,iBeg:iEnd,cp]))
            ydata=cropYield[icity,iBeg:iEnd,cp]
            ydata=np.ma.compressed(np.ma.masked_array(ydata,badYearYield[0,iBeg:iEnd,cp]))
            if size(ydata)==0: # don't plot stuff with no data
                continue
                
            if size(x)==0: # don't plot stuff with no data
                continue
                
            ydataAvg=Avg(ydata)
            slope[iplot,icity,cp],bIntercept[iplot,icity,cp]=polyfit(x,ydata,1)
            yfit=slope[iplot,icity,cp]*x+bIntercept[iplot,icity,cp]
            
            if makePlots:
                figure(1,figsize=(9,4))
                plot(x,ydata,'--*b',x,yfit,'g')
                ylabel('Yield, Bushels/Acre')
                xlabel('year')
                title(ncrop[cp]+' Yield: '+city+' County, IL, slope='+str(round(slope[iplot,icity,cp],2))+' Bu/Acre/Year')
                grid(True)
                savefig('figures/mostcity/'+str(cIndexState[icity])+'/'+city+'/'+ncrop[cp]+'_yield_over_time'+'_little',dpi=700 )
                show()
                clf()
    #0
    iplot+=1 #1
    
    if Normalizedtime: # iplot=1
        
        for cp in range(nCrop):
            if runCrops[cp]==0:
                continue
            x=yr[iBeg:iEnd]
            ydata=cropYield[icity,iBeg:iEnd,cp]
            if size(ydata[badYearYieldBoolean[0,iBeg:iEnd,cp]])==0: # don't plot stuff with no data
                continue
                
            if size(x)==0: # don't plot stuff with no data
                continue
                
            ydataAvg=Avg(ydata[badYearYieldBoolean[0,iBeg:iEnd,cp]])
            slope[iplot,icity,cp],b=polyfit(x[badYearYieldBoolean[0,iBeg:iEnd,cp]],ydata[badYearYieldBoolean[0,iBeg:iEnd,cp]],1)
            yfit=slope[iplot,icity,cp]*x+b
            
            num=len(yfit)-1
            
            cropYield[icity,iBeg:iEnd,cp]=ydata-(slope[iplot,icity,cp]*x+b)
            cropYield[icity,iBeg:iEnd,cp]=cropYield[icity,iBeg:iEnd,cp]+yfit[num]
            
            dataAt2015[icity,cp]=yfit[num]
            
            ydata=np.ma.compressed(np.ma.masked_array(cropYield[icity,iBeg:iEnd,cp],badYearYield[0,iBeg:iEnd,cp]))
            x=np.ma.compressed(np.ma.masked_array(yr[iBeg:iEnd],badYearYield[0,iBeg:iEnd,cp]))
            ydataAvg=Avg(ydata)
            slope[iplot,icity,cp],bIntercept[iplot,icity,cp]=polyfit(x,ydata,1)
            yfit=slope[iplot,icity,cp]*x+bIntercept[iplot,icity,cp]
            Corr[iplot,icity,cp]=corr(x,ydata) 
            
            if makePlots:
                figure(1,figsize=(9,4))
                plot(x,ydata,'*b',x,yfit,'g')
                ylabel('Yield')
                xlabel('year')
                title(city+' '+ncrop[cp]+' Normalized Yield over time  m='+str(round(slope[iplot,icity,cp],3)*100)+
                    ' Corr='+str(round(Corr[iplot,icity,cp],2)))
                grid(True)
                savefig('figures/mostcity/'+str(cIndexState[icity])+'/'+city+'/'+ncrop[cp]+'_normalized_yield_over_time_little',dpi=700)
                show()
                clf()
    #1
    iplot+=1 #2
    
    if UseKDD: # iplot=2
        for cp in range(nCrop):
            if runCrops[cp]==0:
                continue
            x=np.ma.compressed(np.ma.masked_array(KDD[iBeg:iEnd,cp],badYearYield[0,iBeg:iEnd,cp]))
            ydata=cropYield[icity,iBeg:iEnd,cp]
            ydata=ydata[badYearYieldBoolean[0,iBeg:iEnd,cp]]
            if size(ydata)==0: # don't plot stuff with no data
                continue
                
            if size(x)==0: # don't plot stuff with no data
                continue
                   
            ydataAvg=Avg(ydata)
            slope[iplot,icity,cp],bIntercept[iplot,icity,cp]=polyfit(x,ydata,1)
            yfit=slope[iplot,icity,cp]*x+bIntercept[iplot,icity,cp]
            Corr[iplot,icity,cp]=corr(x,ydata)
            R2[iplot,icity,cp]=1.-SumOfSquares(yfit-ydata)/SumOfSquares(ydata-ydataAvg) 
            
            if makePlots:
                figure(1,figsize=(9,4))
                plot(x,ydata,'*b',x,yfit,'g')
                ylabel(ncrop[cp]+' Yield, Bushels/Acre')
                xlabel('Killing Degree Days per Growing Season')
                title(ncrop[cp]+' Yield vs Killing Degree Days, '+city+' IL, Corr='+str(round(Corr[iplot,icity,cp],2)))
                grid(True)
                savefig('figures/mostcity/'+str(cIndexState[icity])+'/'+city+'/'+ncrop[cp]+'_KDD_yield_corr'+'_little',
                    dpi=700 )
                show()
                clf()
    #2
    iplot+=1 #3
    
    if Summer: # iplot=3
        for cp in range(nCrop):
            if runCrops[cp]==0:
                continue
            ydata=cropYield[icity,iBeg:iEnd,cp]
            ydata=np.ma.compressed(np.ma.masked_array(ydata,badYearsforAvg[0,iBeg:iEnd,cp])) 
            x=sAvg[0,iBeg:iEnd,1]
            x=np.ma.compressed(np.ma.masked_array(x,badYearsforAvg[0,iBeg:iEnd,cp]))
            
            if size(ydata)==0: # don't plot stuff with no data
                continue
                
            if size(x)==0: # don't plot stuff with no data
                continue
                
            ydataAvg=Avg(ydata)
            slope[iplot,icity,cp],bIntercept[iplot,icity,cp]=polyfit(x,ydata,1)
            yfit=slope[iplot,icity,cp]*x+bIntercept[iplot,icity,cp]
            Corr[iplot,icity,cp]=corr(x,ydata) 
            R2[iplot,icity,cp]=1.-SumOfSquares(yfit-ydata)/SumOfSquares(ydata-ydataAvg)

            if makePlots:
                figure(1,figsize=(9,4))
                plot(x,ydata,'*b',x,yfit,'g')
                ylabel(ncrop[cp]+' Yield, Bushels/Acre')
                xlabel('Avg Summer (JJA) Temp, F')
                title(ncrop[cp]+' Yield vs Summer Avg Temp, '+city+' IL, Corr='+str(round(Corr[iplot,icity,cp],2)))
                grid(True)
                savefig('figures/mostcity/'+str(cIndexState[icity])+'/'+city+'/'+ncrop[cp]+'_summertemp_yield_corr_'+city+'_little',
                    dpi=700 )
                show()
                clf()
    #3
    iplot+=1 #4
    
              
    if UseHeatWaves: # iplot=4
        for cp in range(nCrop):
            if runCrops[cp]==0:
                continue
            ydata=cropYield[icity,iBeg:iEnd,cp]
            ydata=np.ma.compressed(np.ma.masked_array(ydata,badYearsforAvg[0,iBeg:iEnd,cp]))
            x=seasonheatwavecount[0,iBeg:iEnd,cp]
            x=np.ma.compressed(np.ma.masked_array(x,badYearsforAvg[0,iBeg:iEnd,cp]))
            
            if size(ydata)==0: # don't plot stuff with no data
                continue
                
            if size(x)==0: # don't plot stuff with no data
                continue
                
            ydataAvg=Avg(ydata)
            slope[iplot,icity,cp],bIntercept[iplot,icity,cp]=polyfit(x,ydata,1)
            yfit=slope[iplot,icity,cp]*x+bIntercept[iplot,icity,cp]
            Corr[iplot,icity,cp]=corr(x,ydata) 
            R2[iplot,icity,cp]=1.-SumOfSquares(yfit-ydata)/SumOfSquares(ydata-ydataAvg)
            
            if makePlots:
                figure(1,figsize=(9,4))
                plot(x,ydata,'*b',x,yfit,'g')
                ylabel(ncrop[cp]+' Yield, Bushels/Acre')
                xlabel('Heat Waves, Number/Growing Season')
                title(ncrop[cp]+' Yield vs Heat Waves, '+city+' IL, Corr='+str(round(Corr[iplot,icity,cp],2)))
                grid(True)
                savefig('figures/mostcity/'+str(cIndexState[icity])+'/'+city+'/'+ncrop[cp]+'_heat_waves_yield_corr_'+city+'_little',
                    dpi=700 )
                show()
                clf()

    for cp in range(nCrop):
        if runCrops[cp]==0:
                continue
        for iplot in range(5):
            if iplot==0:
                maxCorr[icity,cp]=Corr[iplot,icity,cp]
                maxCorrAbs=abs(Corr[iplot,icity,cp])
                iplotCorr[icity,cp]=iplot
                maxCorrAbs2=abs(Corr[iplot,icity,cp])
                
            if abs(Corr[iplot,icity,cp])>maxCorrAbs:
                maxCorr2[icity,cp]=maxCorr[icity,cp]
                maxCorrAbs2=abs(maxCorr2[icity,cp])
                iplotCorr2[icity,cp]=iplotCorr[icity,cp]
                
                maxCorr[icity,cp]=Corr[iplot,icity,cp]
                maxCorrAbs=abs(Corr[iplot,icity,cp])
                iplotCorr[icity,cp]=iplot
                
            if maxCorrAbs2<abs(Corr[iplot,icity,cp])<maxCorrAbs:
                maxCorr2[icity,cp]=Corr[iplot,icity,cp]
                maxCorrAbs2=abs(Corr[iplot,icity,cp])
                iplotCorr2[icity,cp]=iplot
                
    if makePlots:
        figure(1,figsize=(9,4))
        #                                   ##
        ## Summer Avg Correlation with Yield ##
        ##                                   ##
        ticks=([110,120,130,140,150,160,170,180,190,200,210,220],[35,40,45,50,55,60])
        
        x=np.arange(1970,2015)
        ydata=sAvg[0,iBeg:iEnd,1]
        
        x2=range(1970,2015)
        ydata2=cropYield[icity,iBeg:iEnd,0]
        
        fig, ax1 = plt.subplots()
        
        ax2 = ax1.twinx()
        ax1.plot(x,ydata,'-*g')
        ax2.plot(x2,ydata2,'-*b')
        ax1.set_yticks([70,75,80,85,90])
        ax2.set_yticks([100,120,140,160,180,200,220,240,260,280])
        ax2.set_ylabel('Detrended Corn Yield, Bushels/Acre',color='b')
        ax1.set_ylabel('Summer Avg Temperature, F',color='g')
        #ax1.set(color='g')
        #ax2.set(color='b')
        ax1.set_xlabel('year')
        title(ncrop[0]+' Yield and Summer Avg Temp, '+city+', IL, Corr='+str(round(Corr[3,icity,0],2)))
        plt.minorticks_on()
        ax1.add_line(Line2D([1980,1980],[50,100],linewidth=.5,color='k',linestyle='dashed'))
        ax1.add_line(Line2D([1983,1983],[50,100],linewidth=.5,color='k',linestyle='dashed'))
        ax1.add_line(Line2D([1988,1988],[50,100],linewidth=.5,color='k',linestyle='dashed'))
        ax1.add_line(Line2D([1991,1991],[50,100],linewidth=.5,color='k',linestyle='dashed'))
        ax1.add_line(Line2D([1995,1995],[50,100],linewidth=.5,color='k',linestyle='dashed'))
        ax1.add_line(Line2D([2002,2002],[50,100],linewidth=.5,color='k',linestyle='dashed'))
        ax1.add_line(Line2D([2012,2012],[50,100],linewidth=.5,color='k',linestyle='dashed'))
        savefig('figures/mostcity/'+str(cIndexState[icity])+'/'+city+'/'+ncrop[0]+'_corr_with_summer_avg_little',dpi=700 )
        show()
        clf()
    
    print 'Corn:'      
    print round(maxCorr[icity,0],2),round(R2[iplotCorr[icity,0],icity,0],2),iplotCorr[icity,0]
    print round(maxCorr2[icity,0],2),round(R2[iplotCorr2[icity,0],icity,0],2),iplotCorr2[icity,0],'\n' 
    print 'Soybeans:'      
    print round(maxCorr[icity,1],2),round(R2[iplotCorr[icity,0],icity,1],2),iplotCorr[icity,1]
    print round(maxCorr2[icity,1],2),round(R2[iplotCorr2[icity,0],icity,1],2),iplotCorr2[icity,1],'\n'
    print 'Rice:'      
    print round(maxCorr[icity,2],2),round(R2[iplotCorr[icity,0],icity,2],2),iplotCorr[icity,2]
    print round(maxCorr2[icity,2],2),round(R2[iplotCorr2[icity,0],icity,2],2),iplotCorr2[icity,2],'\n'

####################################################################
# Write slope, corr, lat, lon, station name into a text document
####################################################################
if not os.path.exists(r'final_data/'):
    os.makedirs(r'final_data/')
pickle.dump(dataAt2015,open('pickle_files/dataAt2015.p','wb'))
pickle.dump(slope,open('pickle_files/slope.p','wb'))
pickle.dump(bIntercept,open('pickle_files/bIntercept.p','wb'))
pickle.dump(lat,open('pickle_files/lat.p','wb'))
pickle.dump(lon,open('pickle_files/lon.p','wb'))
pickle.dump(stationList,open('pickle_files/station.p','wb'))
pickle.dump(goodCity,open('pickle_files/goodCity.p','wb'))
pickle.dump(cityList,open('pickle_files/cityList.p','wb'))
pickle.dump(Corr,open('pickle_files/Corr.p','wb'))
pickle.dump(maxCorr,open('pickle_files/maxCorr.p','wb'))
pickle.dump(iplotCorr,open('pickle_files/iplotCorr.p','wb'))
pickle.dump(maxCorr2,open('pickle_files/maxCorr2.p','wb'))
pickle.dump(iplotCorr2,open('pickle_files/iplotCorr2.p','wb'))
pickle.dump(T90all,open('pickle_files/T90all.p','wb'))
pickle.dump(T10all,open('pickle_files/T10all.p','wb'))
pickle.dump(SummerAvg,open('pickle_files/SummerAvg.p','wb'))
pickle.dump(HeatWaves,open('pickle_files/HeatWaves.p','wb'))
pickle.dump(KDDays,open('pickle_files/KDD.p','wb'))
pickle.dump(cropYield,open('pickle_files/NormalizedCropYield.p','wb'))

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################

###########################################################
#Plots future crop yields
###########################################################
from pylab import *
import csv
from math import sqrt
from sys import exit
import pickle
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from mpl_toolkits.basemap import Basemap
import os
from scipy.stats import norm
import matplotlib as mpl
from matplotlib.patches import Polygon

###############################################
# Functions
###############################################
def Avg(x):   
    '''function to average'''
    xAvg=0.
    for k in range(len(x)):
        xAvg=xAvg+x[k]
    xAvg=xAvg/(k+1)
    return xAvg 
    
    
### Running mean/Moving average
def movingaverage(interval, window_size):
    window = np.ones(int(window_size))/float(window_size)
    return np.convolve(interval, window, 'same')
       
###############################################
# Variables
###############################################
futureYield=pickle.load(open('pickle_files/futureYield.p','rb'))
# futureYield dimensions = (nCities,nPredictor,nyears,nScen,nCrop)
# nPredictor  0=KDD, 1=Summer Avg, 2=Heat Waves
cropYield=pickle.load(open('pickle_files/cropYield.p','rb'))
# cropYield dimensions = (nCities,year,cp)
NormalizedCropYield=pickle.load(open('pickle_files/NormalizedCropYield.p','rb'))
cIndex=pickle.load(open('pickle_files/cIndex.p','rb'))
cIndexState=pickle.load(open('pickle_files/cIndexState.p','rb'))
presentGrowingCounties=pickle.load(open('pickle_files/presentGrowingCounties.p','rb'))
SummerAvg=pickle.load(open('pickle_files/SummerAvg.p','rb'))
HeatWaves=pickle.load(open('pickle_files/HeatWaves.p','rb'))
KDDays=pickle.load(open('pickle_files/KDD.p','rb'))
highTotalYield=pickle.load(open('pickle_files/highTotalYield.p','rb'))
dataAt2015=pickle.load(open('pickle_files/dataAt2015.p','rb'))
slope=pickle.load(open('pickle_files/slope.p','rb'))
bIntercept=pickle.load(open('pickle_files/bIntercept.p','rb'))
SummerAvgFuture=pickle.load(open('pickle_files/SummerAvgFuture.p','rb'))
BBcounty2=pickle.load(open('pickle_files/BBcounty2.p','rb'))

nScen=2
nCrop=3
goodfuture=zeros(shape=(nScen,3143,nCrop),dtype=bool)
goodpast=zeros(shape=(3143,nCrop),dtype=bool)
goodpastyears=zeros(shape=(3143,116,nCrop),dtype=bool)
nyears=100
avgFutureYield=zeros(shape=(nScen,nyears,nCrop))
avgFutureYieldIndividualCities=zeros(shape=(nScen,3143,nyears,nCrop))
avgCropYield=zeros(shape=(116,nCrop))
futureYieldPlusTrend=zeros(shape=(nScen,3143,100,nCrop))
avgFutureYieldPlusTrend=zeros(shape=(nScen,100,nCrop))
avgNormalizedCropYield=zeros(shape=(116,nCrop))  
SummerAvgAvg=zeros(shape=(116,nCrop)) 
goodSummer=zeros(shape=(3143,nCrop),dtype=bool)
goodSummeryears=zeros(shape=(3143,116,nCrop),dtype=bool)

SummerAvgAvg2=zeros(shape=(116)) 
goodSummer2=zeros(shape=(3143),dtype=bool)
goodSummeryears2=zeros(shape=(3143,116),dtype=bool)

SummerAvgAvgFuture=zeros(shape=(2,116)) 
goodSummerFuture=zeros(shape=(2,3143),dtype=bool)
goodSummeryearsFuture=zeros(shape=(2,3143,116),dtype=bool)

cropTitle=('Corn','Soybean','Rice') 
subtract=(20,5,200)
units=('Bushels/Acre','Bushels/Acre','Pounds/Acre')

for scen in range(nScen):
    for y in range(16,100):
        for cp in range(nCrop):
            for icity in range(3143):
                if presentGrowingCounties[icity,cp]==False:
                    continue
                avgFutureYieldIndividualCities[scen,icity,y,cp]=mean(futureYield[icity,0:2,y,scen,cp])
            
yr=arange(2016,2100)
x=yr
for icity in range(3143): 
    if sum(presentGrowingCounties[icity,:])==0 or icity==24:
        continue           
    for cp in range(nCrop):
        if presentGrowingCounties[icity,cp]==False:
            continue
        for scen in range(nScen):
            futureYieldPlusTrend[scen,icity,16:100,cp]=avgFutureYieldIndividualCities[scen,icity,16:100,cp]-dataAt2015[icity,cp]
            futureYieldPlusTrend[scen,icity,16:100,cp]=futureYieldPlusTrend[scen,icity,16:100,cp]+(slope[0,icity,cp]*x+bIntercept[0,icity,cp])+subtract[cp]
        
        
for icity in range(3143):
    for scen in range(nScen):
        for y in range(16,100):
            for cp in range(nCrop):
                if futureYield[icity,0,y,scen,cp]>0 and presentGrowingCounties[icity,cp]==1 and isnan(futureYield[icity,0,y,scen,cp])==False:
                    goodfuture[scen,icity,cp]=True
            if SummerAvgFuture[scen,icity,y]>0 and BBcounty2[icity]==True:
                    goodSummeryearsFuture[scen,icity,y]=True

    for y in range(70,116):
        for cp in range(nCrop):
            if cropYield[icity,y,cp]>0 and highTotalYield[icity,cp]==True:
                    goodpast[icity,cp]=True
                    goodpastyears[icity,y,cp]=True
            if SummerAvg[icity,y]>0 and highTotalYield[icity,cp]==True:
                    goodSummeryears[icity,y,cp]=True
            if SummerAvg[icity,y]>0 and BBcounty2[icity]==True:
                    goodSummeryears2[icity,y]=True
                            
for scen in range(nScen):
    for y in range(16,100):
        for cp in range(nCrop):
            avgFutureYield[scen,y,cp]=mean(futureYield[goodfuture[scen,:,cp],0:2,y,scen,cp])+subtract[cp]
            avgFutureYieldPlusTrend[scen,y,cp]=mean(futureYieldPlusTrend[scen,goodfuture[scen,:,cp],y,cp])
    
j=zeros(shape=(116,nCrop))
summerj=zeros(shape=(116,nCrop))
summerjFuture=zeros(shape=(2,116))
summerj2=zeros(shape=(116))
for cp in range(nCrop):
    for icity in range(3143):
        if goodpast[icity,cp]==False:
            continue
        for y in range(70,116):
            if goodpastyears[icity,y,cp]==True:
                avgCropYield[y,cp]+=cropYield[icity,y,cp]
                avgNormalizedCropYield[y,cp]+=NormalizedCropYield[icity,y,cp]
                j[y,cp]+=1
            if goodSummeryears[icity,y,cp]==True:
                SummerAvgAvg[y,cp]+=SummerAvg[icity,y]
                summerj[y,cp]+=1
                
            if goodSummeryears2[icity,y]==True:
                SummerAvgAvg2[y]+=SummerAvg[icity,y]
                summerj2[y]+=1 
                
        for y in range(16,100):
            for scen in range(2):
                if goodSummeryearsFuture[scen,icity,y]==False:
                    continue
                SummerAvgAvgFuture[scen,y]+=SummerAvgFuture[scen,icity,y]
                summerjFuture[scen,y]+=1                

for y in range(70,116):
    for cp in range(nCrop):
        avgCropYield[y,cp]=avgCropYield[y,cp]/j[y,cp]
        SummerAvgAvg[y,cp]=SummerAvgAvg[y,cp]/summerj[y,cp]
        avgNormalizedCropYield[y,cp]=(avgNormalizedCropYield[y,cp]/j[y,cp])#-subtract[cp]
    SummerAvgAvg2[y]=SummerAvgAvg2[y]/summerj2[y]
for y in range(16,100):
    for scen in range(2):
        SummerAvgAvgFuture[scen,y]=SummerAvgAvgFuture[scen,y]/summerjFuture[scen,y]

###############################################
# Future and Past Summer Avg temp
###############################################
x=range(2016,2100)
ydata=SummerAvgAvgFuture[0,16:100]
yfit=movingaverage(ydata,5)

x2=range(1970,2016)
ydata2=SummerAvgAvg2[70:116]
yfit2=movingaverage(ydata2,5)

x3=range(2016,2100)
ydata3=SummerAvgAvgFuture[1,16:100]
yfit3=movingaverage(ydata3,5)

#plot(x,ydata,'*b',x,yfit,'g',x2,ydata2,'*k',x2,yfit2,'g')
plot(x2[2:44],yfit2[2:44],'-g',linewidth=3)
plot(x3[2:82],yfit3[2:82],'-b',linewidth=3)
plot(x[2:82],yfit[2:82],'-r',linewidth=3)
plot(x2,ydata2,'-g',linewidth=1)
plot(x,ydata,'-r',linewidth=1)
plot(x3,ydata3,'-b',linewidth=1)
legend(['past: historical','future: low emissions','future: high emissions'],loc='lower center')
ylabel(cropTitle[cp]+' Yield, '+units[cp])
xlabel('year')
title('Bread Basket Summer Average Temp: Past vs Future')
grid(True)
savefig('final_figures/poster_final/past_vs_future_summer',dpi=700)
show()

clf()    
    

###############################################
# Past Yields
###############################################

for cp in range(nCrop):
    
    x2=np.arange(1970,2016)
    ydata2=avgCropYield[70:116,cp]
    ydataAvg2=Avg(ydata2)
    slope2,b2=polyfit(x2,ydata2,1)
    yfit2=slope2*x2+b2
    
    plot(x2,ydata2,'--*b',x2,yfit2,'-g')
    ylabel(cropTitle[cp]+' Yield, Bushels/Acre')
    xlabel('year')
    title(cropTitle[cp]+' Yield: Avg of Corn Growing Counties, slope='+str(round(slope2,2))+' Bu/Acre/Year')
    grid(True)
    savefig('final_figures/yield/'+cropTitle[cp]+'/past_yield_'+cropTitle[cp],dpi=700)
    show()
    exit()
    clf()
  
###############################################
# Summer Avg correlation with Yield
###############################################
ticks=([110,120,130,140,150,160,170,180,190,200,210,220],[35,40,45,50,55,60])

for cp in range(nCrop):
    if cp==2:
        continue
    x=np.arange(1970,2016)
    ydata=SummerAvgAvg[70:116,cp]
    #ydataAvg=Avg(ydata)
    #slope,b=polyfit(x,ydata,1)
    #yfit=slope*x+b
    
    x2=range(1970,2016)
    ydata2=avgNormalizedCropYield[70:116,cp]
    #ydataAvg2=Avg(ydata2)
    #slope2,b2=polyfit(x2,ydata2,1)
    #yfit2=slope2*x2+b2
    
    fig, ax1 = plt.subplots()
    
    ax2 = ax1.twinx()
    ax1.plot(x,ydata,'-*g')
    ax2.plot(x2,ydata2,'-*b')
    ax1.set_yticks([70,75,80,85,90])
    ax2.set_yticks(ticks[cp])
    ax2.set_ylabel(cropTitle[cp]+' Yield',color='b')
    ax1.set_ylabel('Summer Avg Temperature',color='g')
    ax1.set_xlabel('year')
    title(cropTitle[cp]+' Yield and Summer Avg Temperature')
    savefig('final_figures/yield/'+cropTitle[cp]+'/'+cropTitle[cp]+'yield_and_summer_avg_temp',dpi=700)
    show()
    clf()    


###############################################
# Future and Past Yields
###############################################
for cp in range(3):
    #figure(1,figsize=(9,5))
    x=np.arange(2016,2100)
    ydata=avgFutureYield[0,16:100,cp]
    yfit=movingaverage(ydata,5)
    slope,bIntercept=polyfit(x,ydata,1)
    yfitline=slope*x+bIntercept
    
    x2=np.arange(1970,2016)
    ydata2=avgCropYield[70:116,cp]
    yfit2=movingaverage(ydata2,5)
    slope2,bIntercept2=polyfit(x2,ydata2,1)
    yfitline2=slope2*x2+bIntercept2
    
    x3=np.arange(2016,2100)
    ydata3=avgFutureYield[1,16:100,cp]
    yfit3=movingaverage(ydata3,5)
    slope3,bIntercept3=polyfit(x3,ydata3,1)
    yfitline3=slope3*x3+bIntercept3
    
    #plot(x,ydata,'*b',x,yfit,'g',x2,ydata2,'*k',x2,yfit2,'g')
    plot(x2[2:44],yfit2[2:44],'-g',linewidth=3)
    plot(x3[2:82],yfit3[2:82],'-b',linewidth=3)
    plot(x[2:82],yfit[2:82],'-r',linewidth=3)
    
    plot(x2,ydata2,'-g',linewidth=1)
    plot(x,ydata,'-r',linewidth=1)
    plot(x3,ydata3,'-b',linewidth=1)
    
    plot(x2,yfitline2, '--g',linewidth=2)
    plot(x,yfitline, '--r',linewidth=2)
    plot(x3,yfitline3, '--b',linewidth=2)
    legend(['past: historical','future: low emissions','future: high emissions'],loc='lower center')
    ylabel(cropTitle[cp]+' Yield, '+units[cp])
    xlabel('year')
    title('U.S. '+cropTitle[cp]+' Yield: Past vs Future')
    grid(True)
    #savefig('final_figures/yield/'+cropTitle[cp]+'/past_vs_future_yield_'+cropTitle[cp]+'_little',dpi=700)
    show()
    print cropTitle[cp],'\nhistorical:',round((slope2*10/yfitline2[0])*100,2),'\nfuture high:',round((slope*10/yfitline[0])*100,2), '\nfuture low:',round((slope3*10/yfitline3[0])*100,2)
    
    pause(1)
    clf() 

###############################################
# Future and Normalized Past Yields
###############################################
for cp in range(nCrop):
    
    x=range(2016,2100)
    ydata=avgFutureYield[0,16:100,cp]
    yfit=movingaverage(ydata,5)
    
    x2=range(1970,2016)
    ydata2=avgNormalizedCropYield[70:116,cp]
    yfit2=movingaverage(ydata2,5)
    
    x3=range(2016,2100)
    ydata3=avgFutureYield[1,16:100,cp]
    yfit3=movingaverage(ydata3,5)
    
    #plot(x,ydata,'*b',x,yfit,'g',x2,ydata2,'*k',x2,yfit2,'g')
    plot(x2[2:44],yfit2[2:44],'-g',linewidth=3)
    plot(x3[2:82],yfit3[2:82],'-b',linewidth=3)
    plot(x[2:82],yfit[2:82],'-r',linewidth=3)
    plot(x2,ydata2,'-g',linewidth=1)
    plot(x,ydata,'-r',linewidth=1)
    plot(x3,ydata3,'-b',linewidth=1)
    legend(['past: detrended','future: low emissions','future: high emissions'],loc='lower center')
 
    ylabel(cropTitle[cp]+' Yield, '+units[cp])
    xlabel('year')
    title('U.S. '+cropTitle[cp]+' Yield: Past vs Future')
    grid(True)
    savefig('final_figures/yield/'+cropTitle[cp]+'/normalized_past_vs_future_yield_'+cropTitle[cp],dpi=700)
    show()
    pause(5)
    clf()    

###############################################
# Future and Past Yields with added in trend
###############################################

for cp in range(nCrop):
    
    x=range(2016,2100)
    ydata=avgFutureYieldPlusTrend[0,16:100,cp]
    yfit=movingaverage(ydata,5)
    
    
    x4=np.arange(2018,2100)
    ydata2=avgCropYield[70:116,cp]
    yfit2=movingaverage(ydata2,5)
    x2=np.arange(1970,2016)
    slope,bIntercept=polyfit(x2,ydata2,1)
    yfit21=slope*x4+bIntercept
    
    x3=range(2016,2100)
    ydata3=avgFutureYieldPlusTrend[1,16:100,cp]
    yfit3=movingaverage(ydata3,5)
    
    #plot(x,ydata,'*b',x,yfit,'g',x2,ydata2,'*k',x2,yfit2,'g')
    plot(x2[2:44],yfit2[2:44],'-g',linewidth=3)

    plot(x3[2:82],yfit3[2:82],'-b',linewidth=3)
    plot(x[2:82],yfit[2:82],'-r',linewidth=3)
    plot(x2,ydata2,'-g',linewidth=1)
    plot(x,ydata,'-r',linewidth=1)
    plot(x3,ydata3,'-b',linewidth=1)
    plot(x4,yfit21,'--g',linewidth=3)
    legend(['past: historical','future: low emissions with trend','future: high emissions with trend'],loc='lower right')
    ylabel(cropTitle[cp]+' Yield, '+units[cp])
    xlabel('year')
    title('U.S. '+cropTitle[cp]+' Yield: Past vs Future')
    grid(True)
    savefig('final_figures/yield/'+cropTitle[cp]+'/past_vs_future_yield_plus_trend_'+cropTitle[cp]+'_little',dpi=700)
    show()
    pause(3)
    clf()

###############################################
# Past Yields US map
###############################################
for cp in range(nCrop):
    #cmapArray=concatenate((plt.cm.jet(arange(128)),plt.cm.PuOr(arange(128))),axis=0)
    j=-1
    MinMaxArray=ones(shape=(2,1))
    
    cmapArray=plt.cm.jet(arange(256))
    cmin=0
    y1=0
    y2=256
    cmax=.8*(np.amax(cropYield[:,105:116,cp]))
    
    subPlot1 = plt.axes([.5,.2,.475,.6])
    MinMaxArray[0,0]=cmin
    MinMaxArray[1,0]=cmax
    plt.imshow(MinMaxArray,cmap='jet')
    plt.colorbar()
    
    # create the map
    subPlot1 = plt.axes([0.1,0,0.8,1])
    map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
    
    # load the shapefile, use the name 'states'
    map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
    
    ax = plt.gca() # get current axes instance
    
    for shape_dict in map.states_info:
        j+=1
        seg = map.states[j]
        
        s=int(shape_dict['STATEFP'])
        c=int(shape_dict['COUNTYFP'])
        
        if s==72 or cIndex[s,c]==-9999:
            continue
    
        Yield=mean(cropYield[cIndex[s,c],70:81,cp])
        
        if Yield<=0:
            continue
            
        x=Yield
        y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
        icmap=min(255,int(round(y,1)))
        
        poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
        ax.add_patch(poly)
    
    title('Avg '+units[cp]+' of '+cropTitle[cp]+' by County 1970-1980')       
    savefig('final_figures/yield/'+cropTitle[cp]+'/'+'1970-1980_yield_us_map_'+cropTitle[cp],dpi=700)
    show()
    clf()

###############################################
# Present Yields US map
###############################################
for cp in range(nCrop):
    j=-1
    MinMaxArray=ones(shape=(2,1))
    
    cmapArray=plt.cm.jet(arange(256))
    cmin=0
    y1=0
    y2=256
    cmax=.8*(np.amax(cropYield[:,105:116,cp]))
    
    subPlot1 = plt.axes([.5,.2,.475,.6])
    MinMaxArray[0,0]=cmin
    MinMaxArray[1,0]=cmax
    plt.imshow(MinMaxArray,cmap='jet')
    plt.colorbar()
    
    # create the map
    subPlot1 = plt.axes([0.1,0,0.8,1])
    map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
    
    # load the shapefile, use the name 'states'
    map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
    
    ax = plt.gca() # get current axes instance
    
    for shape_dict in map.states_info:
        j+=1
        seg = map.states[j]
        
        s=int(shape_dict['STATEFP'])
        c=int(shape_dict['COUNTYFP'])
        
        if s==72 or cIndex[s,c]==-9999:
            continue
        
        Yield=mean(cropYield[cIndex[s,c],105:116,cp])
        
        if Yield<=0:
            continue
            
        x=Yield
        y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
        icmap=min(255,int(round(y,1)))
        
        poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
        ax.add_patch(poly)
    
    title('Avg '+units[cp]+' of '+cropTitle[cp]+' by County 2005-2015')       
    savefig('final_figures/yield/'+cropTitle[cp]+'/'+'2005-2015_yield_us_map_'+cropTitle[cp],dpi=700)
    show() 
    clf()
 
###############################################
# Future Yields US map
###############################################
scen=0
for cp in range(nCrop):
    j=-1 #counter for how many times loop goes through
    
    MinMaxArray=ones(shape=(2,1))
    cmapArray=plt.cm.jet(arange(256))
    cmin=0
    y1=0
    y2=256
    cmax=.8*(np.amax(cropYield[:,105:116,cp]))
    
    subPlot1 = plt.axes([.5,.2,.475,.6])
    MinMaxArray[0,0]=cmin
    MinMaxArray[1,0]=cmax
    plt.imshow(MinMaxArray,cmap='jet')
    plt.colorbar()
    
    # create the map
    subPlot1 = plt.axes([0.1,0,0.8,1])
    map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
    
    # load the shapefile, use the name 'states'
    map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)
    
    ax = plt.gca() # get current axes instance
    
    for shape_dict in map.states_info:
        j+=1
        seg = map.states[j]
        
        s=int(shape_dict['STATEFP'])
        c=int(shape_dict['COUNTYFP'])
    
        if s==72 or cIndex[s,c]==-9999:
            continue
        
        Yield=mean(futureYield[cIndex[s,c],:,90:100,scen,cp])    
        
        if isnan(Yield):
            continue
        
        if Yield<=0:
            continue
            
        x=Yield
        y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
        icmap=min(255,int(round(y,1)))
        
        poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
        ax.add_patch(poly)
    
    title('Avg '+units[cp]+' of '+cropTitle[cp]+' by County 2090-2100, High Emissions')       
    savefig('final_figures/yield/'+cropTitle[cp]+'/'+'2090-2100_yield_us_map_'+cropTitle[cp],dpi=700)
    show() 
    clf() 
        
    
    
    
#presentGrowingCounties=zeros(shape=(3143,3))
#for icity in range(3143):
#    for cp in range(3):
#        Yield=mean(cropYield[icity,105:115,cp])
#        if Yield>0 or highTotalYield[icity,cp]==True:
#            presentGrowingCounties[icity,cp]=1

    
###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################    
    
###########################################################
# Plots future crop yields
###########################################################
from pylab import *
import csv
from math import sqrt
from sys import exit
import pickle
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from mpl_toolkits.basemap import Basemap
import os
from scipy.stats import norm
import matplotlib as mpl
from matplotlib.patches import Polygon

###############################################
# Functions
###############################################
def Avg(x):   
    '''function to average'''
    xAvg=0.
    for k in range(len(x)):
        xAvg=xAvg+x[k]
    xAvg=xAvg/(k+1)
    return xAvg 
    
    

futureYield=pickle.load(open('pickle_files/futureYield.p','rb'))
# futureYield dimensions = (nCities,nPredictor,nyears,nScen,nCrop)
# nPredictor  0=Summer avg, 1=Heat Waves, 2=KDD
cropYield=pickle.load(open('pickle_files/cropYield.p','rb'))
# cropYield dimensions = (nCities,year,cp)
cIndex=pickle.load(open('pickle_files/cIndex.p','rb'))
presentGrowingCounties=pickle.load(open('pickle_files/presentGrowingCounties.p','rb'))

nyears=32
cp=0

icity=789
avgFutureYield=zeros(shape=(nyears))
avgCropYield=zeros(shape=(116))        

for y in range(nyears):
    avgFutureYield[y]=mean(futureYield[presentGrowingCounties[cp],0:2,y,0,cp])
for y in range(70,116):
    avgCropYield[y]=mean(cropYield[presentGrowingCounties[cp],y,cp])


x=range(2020,2052)
ydata=avgFutureYield 
#ydataAvg=Avg(ydata)
#slope,b=polyfit(x,ydata,1)
#yfit=slope*x+b

x2=range(1970,2016)
ydata2=avgCropYield[70:116]
#ydataAvg2=Avg(ydata2)
#slope2,b2=polyfit(x2,ydata2,1)
#yfit2=slope2*x2+b2

#plot(x,ydata,'*b',x,yfit,'g',x2,ydata2,'*k',x2,yfit2,'g')
plot(x2,ydata2,'-*g',x,ydata,'-*b')
legend(['past yield','future yield'],loc='upper left')
ylabel('Yield')
xlabel('year')
title('Past vs Future Yield')
grid(True)
show()
clf()

# create the map
map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
        projection='lcc',lat_1=33,lat_2=45,lon_0=-95)

# load the shapefile, use the name 'states'
map.readshapefile('data/shape_files/counties/cb_2015_us_county_20m', name='states', drawbounds=True)

ax = plt.gca() # get current axes instance
j=0
cmapArray=plt.cm.jet(arange(256))
cmin=0
y1=0
y2=256
year=100
cp=2
cmax=np.amax(cropYield[:,year,cp])

for shape_dict in map.states_info:
    seg = map.states[j]
    
    s=int(shape_dict['STATEFP'])
    c=int(shape_dict['COUNTYFP'])
    
    if s==72 or cIndex[s,c]==-9999:
        j+=1
        continue

    Yield=cropYield[cIndex[s,c],year,cp]
    
    if Yield==-9999.0 or Yield==0:
        j+=1
        continue
        
    x=Yield
    y=y1+(y2-y1)/(cmax-cmin)*(x-cmin)
    icmap=min(255,int(round(y,1)))
       
    poly = Polygon(seg,facecolor=[cmapArray[icmap,0],cmapArray[icmap,1],cmapArray[icmap,2]],edgecolor=[0,0,0])
    ax.add_patch(poly)
       
    j+=1
title('Crop '+str(cp)+' Yield for Year '+str(year+1900))
plt.show() 
# the Scientific Python netCDF 3 interface
# http://dirac.cnrs-orleans.fr/ScientificPython/
# from Scientific.IO.NetCDF import NetCDFFile as Dataset
# the 'classic' version of the netCDF4 python interface
# http://code.google.com/p/netcdf4-python/
from netCDF4 import Dataset
from numpy import arange # array module from http://numpy.scipy.org
from numpy.testing import assert_array_equal, assert_array_almost_equal
from pylab import *
import matplotlib
import csv
from math import sqrt
from sys import exit
import numpy as np
import os
import pickle
from datetime import datetime
from datetime import timedelta
import sys

TMAX=True
TMIN=False

rcp85=False
rcp45=True

# open a the netCDF file for reading.
if rcp85:
    if TMAX:    
        ncfile = (Dataset('tasmax_rcp85_2092_2099.nc','r'),
            Dataset('tasmax_rcp85_2084_2091.nc','r'),
            Dataset('tasmax_rcp85_2076_2083.nc','r'),
            Dataset('tasmax_rcp85_2068_2075.nc','r'),
            Dataset('tasmax_rcp85_2060_2067.nc','r'),
            Dataset('tasmax_rcp85_2052_2059.nc','r'),
            Dataset('tasmax_rcp85_2044_2051.nc','r'),
            Dataset('tasmax_rcp85_2036_2043.nc','r'),
            Dataset('tasmax_rcp85_2028_2035.nc','r'),
            Dataset('tasmax_rcp85_2020_2027.nc','r'),
            Dataset('tasmax_rcp85_2016_2019.nc','r'))
        t='tmax'
    
    if TMIN:
        ncfile = (Dataset('tasmin_rcp85_2092_2099.nc','r'),
            Dataset('tasmin_rcp85_2084_2091.nc','r'),
            Dataset('tasmin_rcp85_2076_2083.nc','r'),
            Dataset('tasmin_rcp85_2068_2075.nc','r'),
            Dataset('tasmin_rcp85_2060_2067.nc','r'),
            Dataset('tasmin_rcp85_2052_2059.nc','r'),
            Dataset('tasmin_rcp85_2044_2051.nc','r'),
            Dataset('tasmin_rcp85_2036_2043.nc','r'),
            Dataset('tasmin_rcp85_2028_2035.nc','r'),
            Dataset('tasmin_rcp85_2020_2027.nc','r'),
            Dataset('tasmin_rcp85_2016_2019.nc','r'))
        t='tmin'
        
if rcp45:
    if TMAX:    
        ncfile = (Dataset('tasmax_rcp45_2092_2099.nc','r'),
            Dataset('tasmax_rcp45_2084_2091.nc','r'),
            Dataset('tasmax_rcp45_2076_2083.nc','r'),
            Dataset('tasmax_rcp45_2068_2075.nc','r'),
            Dataset('tasmax_rcp45_2060_2067.nc','r'),
            Dataset('tasmax_rcp45_2052_2059.nc','r'),
            Dataset('tasmax_rcp45_2044_2051.nc','r'),
            Dataset('tasmax_rcp45_2036_2043.nc','r'),
            Dataset('tasmax_rcp45_2028_2035.nc','r'),
            Dataset('tasmax_rcp45_2020_2027.nc','r'),
            Dataset('tasmax_rcp45_2016_2019.nc','r'))
        t='tmax'
    
    if TMIN:
        ncfile = (Dataset('tasmin_rcp45_2092_2099.nc','r'),
            Dataset('tasmin_rcp45_2084_2091.nc','r'),
            Dataset('tasmin_rcp45_2076_2083.nc','r'),
            Dataset('tasmin_rcp45_2068_2075.nc','r'),
            Dataset('tasmin_rcp45_2060_2067.nc','r'),
            Dataset('tasmin_rcp45_2052_2059.nc','r'),
            Dataset('tasmin_rcp45_2044_2051.nc','r'),
            Dataset('tasmin_rcp45_2036_2043.nc','r'),
            Dataset('tasmin_rcp45_2028_2035.nc','r'),
            Dataset('tasmin_rcp45_2020_2027.nc','r'),
            Dataset('tasmin_rcp45_2016_2019.nc','r'))
        t='tmin'
    
# read the data in variable named 'data'.
print 'read lat:'
lat = ncfile[0].variables['lat'][:]
print 'read lon'
lon = ncfile[0].variables['lon'][:]
print 'open nearest station file:'
f=open('../written_files/model_nearest_stations.txt','r')

nCities=3143

modelData45=-9999*ones(shape=(31,12,85,3143),dtype=np.float16)
cmiplat=zeros(shape=(nCities),dtype=np.float16)
cmiplon=zeros(shape=(nCities),dtype=np.float16)

for files in range(11):
    f.close()
    f=open('../written_files/model_nearest_stations.txt','r')
    print '\n\n'+t; sys.stdout.flush()
    print str(files) +' of 10'; sys.stdout.flush()
    time=ncfile[files].variables['time']
            
    icity=-1
    for cityline in f:
        icity+=1 
        print icity,'of',nCities; sys.stdout.flush()
        
        tmp=cityline.split(',')
        lat=tmp[0]
        lon=tmp[1]
        ilat=tmp[2]
        ilon=tmp[3]
        cmiplat[icity]=lat
        cmiplon[icity]=lon
        
        if int(float(lat))==-9999:
            print 'lat='+lat+': continued'; sys.stdout.flush()
            continue   
        if TMAX:
            tmax=ncfile[files].variables['air_temperature'][:,ilat,ilon]
        if TMIN:
            tmin=ncfile[files].variables['air_temperature'][:,ilat,ilon]
        
        for iday in range(len(time)):
            date_format = "%m/%d/%Y"
            d = timedelta(days=int(time[iday]))
            a = datetime.strptime('1/1/1900', date_format)
            date=a+d
            year=date.year
            month=date.month
            day=date.day
            
            y=year-2016
            m=month-1
            d=day-1
            
            if TMAX:    
                modelData45[d,m,y,icity]=(tmax[iday]-273.15)*1.8+32.0 # convert K to F
            if TMIN:
                modelData45[d,m,y,icity]=(tmin[iday]-273.15)*1.8+32.0 # convert K to F 
     
print cmiplat
print cmiplon
print np.shape(cmiplat)

pickle.dump(cmiplat,open('cmiplatall45.p','wb'))
pickle.dump(cmiplon,open('cmiplonall45.p','wb'))
pickle.dump(modelData45,open('modelData_all_'+t+'_rcp45.p','wb'))

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################

###############################################
# Finds growing seasons
###############################################

from pylab import *
import csv
from math import sqrt
from sys import exit
import pickle
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from mpl_toolkits.basemap import Basemap
import os
from scipy.stats import norm
import matplotlib as mpl
from matplotlib.patches import Polygon
import random

# cd Documents/Science_Fair_2017_Crop_Yields/

f=open('data/growing_season_corn.csv','r')

startGrowingMon=-9999*ones(shape=(57,3))
endGrowingMon=-9999*ones(shape=(57,3))
startGrowingDay=-9999*ones(shape=(57,3))
endGrowingDay=-9999*ones(shape=(57,3))

for cp in range(3):
    for character in f:
        tmp=character[0:19]
        
        s=int(tmp[3:5])
        startDate=tmp[6:12]
        endDate=tmp[13:19]
        
        if startDate[0:3]=='Apr':
            startMon=3
        if startDate[0:3]=='May':
            startMon=4
        if startDate[0:3]=='Jun':
            startMon=5
        if startDate[0:3]=='Jul':
            startMon=6
        
        if endDate[0:3]=='Aug':
            endMon=7
        if endDate[0:3]=='Sep':
            endMon=8
        if endDate[0:3]=='Oct':
            endMon=9
        if endDate[0:3]=='Nov':
            endMon=10
        
        startGrowingMon[s,cp]=int(startMon)
        startGrowingDay[s,cp]=int(startDate[4:6])
        endGrowingMon[s,cp]=int(endMon)
        endGrowingDay[s,cp]=int(endDate[4:6])
    if cp==0:
        f=open('data/growing_season_soy.csv','r')
    if cp==1:
        f=open('data/growing_season_rice.csv','r')
    
pickle.dump(startGrowingMon,open('pickle_files/startGrowingMon.p','wb'))
pickle.dump(startGrowingDay,open('pickle_files/startGrowingDay.p','wb'))
pickle.dump(endGrowingMon,open('pickle_files/endGrowingMon.p','wb'))
pickle.dump(endGrowingDay,open('pickle_files/endGrowingDay.p','wb'))

###################################################################################################
########################################## NEW CODE ###############################################
###################################################################################################

#####################################################################################
# Plots different variables for years 1970-80 and 2090-2100: low and high emmisions
#####################################################################################
from pylab import *
import csv
from math import sqrt
from sys import exit
import pickle
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
from mpl_toolkits.basemap import Basemap
import os
from scipy.stats import norm
import matplotlib as mpl

presentGrowingCounties=pickle.load(open('pickle_files/presentGrowingCounties.p','rb'))
BBcounty=pickle.load(open('pickle_files/BBcounty2.p','rb'))
BBsummerPast=pickle.load(open('pickle_files/BBsummerPast.p','rb'))
BBHeatWavesPast=pickle.load(open('pickle_files/BBHeatWavesPast.p','rb'))
BBKDDpast=pickle.load(open('pickle_files/BBKDDpast.p','rb'))
SummerAvgFuture=pickle.load(open('pickle_files/SummerAvgFuture.p','rb'))
KDDFuture=pickle.load(open('pickle_files/KDDFuture.p','rb'))
HeatWavesFuture=pickle.load(open('pickle_files/HeatWavesFuture.p','rb'))

BBsummerFuture85=-9999*ones(shape=(3143,10))
BBHeatWavesFuture85=-9999*ones(shape=(3143,10))
BBKDDfuture85=-9999*ones(shape=(3143,10))
BBsummerFuture45=-9999*ones(shape=(3143,10))
BBHeatWavesFuture45=-9999*ones(shape=(3143,10))
BBKDDfuture45=-9999*ones(shape=(3143,10))
for icity in range(3143):
    if BBcounty[icity]==False:
        continue
    for y in range(10):
        BBsummerFuture85[icity,y]=SummerAvgFuture[0,icity,y+90]
        BBHeatWavesFuture85[icity,y]=HeatWavesFuture[0,icity,y+90]
        BBKDDfuture85[icity,y]=KDDFuture[0,icity,y+90]
        
        BBsummerFuture45[icity,y]=SummerAvgFuture[1,icity,y+90]
        BBHeatWavesFuture45[icity,y]=HeatWavesFuture[1,icity,y+90]
        BBKDDfuture45[icity,y]=KDDFuture[1,icity,y+90]              
       
BBsummerPast=BBsummerPast[:,0:9].ravel()
BBHeatWavesPast=BBHeatWavesPast[:,0:9].ravel()
BBKDDpast=BBKDDpast[:,0:9].ravel()

BBsummerFuture85=BBsummerFuture85.ravel()
BBHeatWavesFuture85=BBHeatWavesFuture85.ravel()
BBKDDfuture85=BBKDDfuture85.ravel()

BBsummerFuture45=BBsummerFuture45.ravel()
BBHeatWavesFuture45=BBHeatWavesFuture45.ravel()
BBKDDfuture45=BBKDDfuture45.ravel()

goodCityFutureLong85Summer=zeros(shape=(62860),dtype=bool)
goodCityFutureLong45Summer=zeros(shape=(62860),dtype=bool)
goodCityPastLongSummer=zeros(shape=(62860),dtype=bool)

goodCityFutureLong85HeatWaves=zeros(shape=(62860),dtype=bool)
goodCityFutureLong45HeatWaves=zeros(shape=(62860),dtype=bool)
goodCityPastLongHeatWaves=zeros(shape=(62860),dtype=bool)

goodCityFutureLong85KDD=zeros(shape=(62860),dtype=bool)
goodCityFutureLong45KDD=zeros(shape=(62860),dtype=bool)
goodCityPastLongKDD=zeros(shape=(62860),dtype=bool)

j=0
for k in range(28287):
    if BBsummerPast[k]>0:
        goodCityPastLongSummer[k]=True
        j+=1
    if BBsummerFuture85[k]>0:
        goodCityFutureLong85Summer[k]=True
    if BBsummerFuture45[k]>0:
        goodCityFutureLong45Summer[k]=True
    
    if BBHeatWavesPast[k]>0:
        goodCityPastLongHeatWaves[k]=True
        j+=1
    if BBHeatWavesFuture85[k]>0:
        goodCityFutureLong85HeatWaves[k]=True
    if BBHeatWavesFuture45[k]>0:
        goodCityFutureLong45HeatWaves[k]=True
    
    if BBKDDpast[k]>0:
        goodCityPastLongKDD[k]=True
        j+=1
    if BBKDDfuture85[k]>0:
        goodCityFutureLong85KDD[k]=True
    if BBKDDfuture45[k]>0:
        goodCityFutureLong45KDD[k]=True

BBvar=(BBsummerPast,BBHeatWavesPast,BBKDDpast,BBsummerFuture85,BBHeatWavesFuture85,BBKDDfuture85,BBsummerFuture45,BBHeatWavesFuture45,BBKDDfuture45)
var=('Summer_Avg_Temp','Heat_Waves','KDD','Summer_Avg_Temp','Heat_Waves','KDD')
labelVar=('Summer Avg Temperature','Heat Waves','Killing Degree Days','Summer Avg Temp','Heat Waves','KDD')
titleVar=('Summer Average Temperature for Bread Basket', #0
    'Heat Waves for Bread Basket',#1
    'Killing Degree Days for Bread Basket') 
binsize=(0.5,1,50,0.5,1,50,0.5,1,50)
minext=(80,0,100,80,0,100,80,0,100)
maxext=(105,50,2500,105,50,2500,105,50,2500)
goodCityLong=(goodCityPastLongSummer,goodCityPastLongHeatWaves,goodCityPastLongKDD,
    goodCityFutureLong85Summer,goodCityFutureLong85HeatWaves,goodCityFutureLong85KDD,
    goodCityFutureLong45Summer,goodCityFutureLong45HeatWaves,goodCityFutureLong45KDD)
ymax=(500,1000,750)

if not os.path.exists(r'final_figures/hists_different_times'):
    os.makedirs(r'final_figures/hists_different_times')
    
for k in range(3):
    figure(1,figsize=(9,5))
    #           #
    # Histogram #
    #           #
    bins1=np.linspace(minext[k],maxext[k],(maxext[k]-minext[k])/binsize[k]+1)
    (mu1, sigma1) = norm.fit(BBvar[k][goodCityLong[k]]) # fit the best fit normal distrobution
    
    bins3=np.linspace(minext[k+6],maxext[k+6],(maxext[k+6]-minext[k+6])/binsize[k+6]+1)
    (mu3, sigma3) = norm.fit(BBvar[k+6][goodCityLong[k+6]]) # fit the best fit normal distrobution
    
    bins2=np.linspace(minext[k+3],maxext[k+3],(maxext[k+3]-minext[k+3])/binsize[k+3]+1)
    (mu2, sigma2) = norm.fit(BBvar[k+3][goodCityLong[k+3]]) # fit the best fit normal distrobution
    
    # plot the histogram
    n, bins, patches = plt.hist(BBvar[k][goodCityLong[k]], bins1, normed=False, histtype='bar', 
        facecolor='green', alpha=0.6, label='1970-1980')
    n, bins, patches = plt.hist(BBvar[k+6][goodCityLong[k+6]], bins3, normed=False, histtype='bar', 
        facecolor='blue', alpha=0.6, label='2090-2100, low emissions')
    n, bins, patches = plt.hist(BBvar[k+3][goodCityLong[k+3]], bins2, normed=False, histtype='bar', 
        facecolor='red', alpha=0.6, label='2090-2100, high emissions')
    plt.legend(loc='upper center')
    if k==0:
        plt.xlabel(labelVar[k])
    else:
        plt.xlabel('Number of '+labelVar[k]+'/year')
    plt.ylabel('Number of instances (years and counties)')
    plt.title(titleVar[k])
    y1 = mlab.normpdf( bins1, mu1, sigma1) # best fit normal ditrobution
    y3 = mlab.normpdf( bins3, mu3, sigma3) # best fit normal ditrobution
    y2 = mlab.normpdf( bins2, mu2, sigma2) # best fit normal ditrobution
    l1 = plt.plot(bins1, y1*binsize[k]*np.size(BBvar[k][goodCityLong[k]]), 'r-', linewidth=2)
    l3 = plt.plot(bins3, y3*binsize[k+6]*np.size(BBvar[k+6][goodCityLong[k+6]]), 'r-', linewidth=2)
    l2 = plt.plot(bins2, y2*binsize[k+3]*np.size(BBvar[k+3][goodCityLong[k+3]]), 'r-', linewidth=2)
    plt.axvline(BBvar[k][goodCityLong[k]].mean(), color='r', linewidth=4) # line at the mean
    plt.axvline(BBvar[k+6][goodCityLong[k+6]].mean(), color='r', linewidth=4) # line at the mean
    plt.axvline(BBvar[k+3][goodCityLong[k+3]].mean(), color='r', linewidth=4) # line at the mean
    axes = plt.gca()
    axes.set_ylim([0,ymax[k]])
    grid(True)
    plt.savefig('final_figures/hists_different_times/'+var[k]+'_overlap_times_little',dpi=500)
    plt.show()
    pause(.5)
    clf()