from pylab import * import csv from math import sqrt #cd /Users/lilli_000/Documents/Supercomputing cityList=['Los_Alamos','Louisiana','Albuquerque_1950','Nebraska', 'California','Cincinnati','Colorado_Springs','Leadville_Colorado', 'New_York_L','North_Dakota','Oregon','Palisade_Colorado', 'ROSWELL','Salt_Lake_City_Utah','Seattle_Washington', 'South_California','South_Carolina'] #startDate=1915 #central valley and New York #startDate=1895 #southern desert startDate=1912 for icity in range(len(cityList)): city=cityList[icity] fileID=open('data/'+city+'.csv','rb') reader=csv.reader(fileID) a=csv.reader j=0 #initalising list variables Tmax=[] Tmin=[] snowd=[] snow=[] precip=[] date=[] year=[] month=[] day=[] #Initialise monthly averages TmaxAvgM=[0.] TminAvgM=[0.] snowdAvgM=[] snowAvgM=[0.] precipAvgM=[0.] dateAvgM=[] yearAvgM=[] monthAvgM=[] #Initialise yearly averages TmaxAvgY=[0.] TminAvgY=[0.] snowdAvgY=[] snowAvgY=[0.] precipAvgY=[0.] dateAvgY=[] yearAvgY=[] monthAvgY=[] #loops over all the rows for row in reader: if j==0: header=row else: tmp=row[1] ytmp=float(tmp[0:4]) #year mtmp=float(tmp[4:6]) #month if ytmp==2015: #choose upper limit year break if ytmp>=startDate: #choose lower limit year #read data into python #divide by 10 to convert 10ths degrees C into degrees C Tmax.append(float(row[5])/10.) #Tmax.append(float(row[4])/10.) Tmin.append(float(row[6])/10.) #Tmin.append(float(row[5])/10.) snow.append(float(row[4])/10./2.54) #snow.append(float(row[3])/10./2.54) precip.append(float(row[2])/100./2.54) tmp=row[1] year.append(float(tmp[0:4])) month.append(float(tmp[4:6])) day.append(float(tmp[6:8])) j=j+1 for k in range(len(year)): date.append(year[k]+(month[k]-1)/12.0+(day[k]-1)/365) #make a date if k>365: interval=365 else: interval=1 if Tmax[k]<-50: #if no data, assume it is yesterday's temperature Tmax[k]=Tmax[k-interval] #print "Tmax = -999" if Tmin[k]<-50: Tmin[k]=Tmin[k-interval] #print "Tmin = -999" if precip[k]<-1: precip[k]=precip[k-interval] #print "precip = -999" if snow[k]<-5: snow[k]=0 #print "snow = -999" #Taking an average of each month dateAvgM.append(year[0]+(month[0]-.5)/12) n=0 #day counter im=0 #month counter for k in range(len(year)): #loop over days if month[k]!=month[k-1] and k>0: #test for new month TmaxAvgM[im]=TmaxAvgM[im]/n TminAvgM[im]=TminAvgM[im]/n n=0 TmaxAvgM.append(0.) TminAvgM.append(0.) precipAvgM.append(0.) snowAvgM.append(0.) dateAvgM.append(year[k]+(month[k]-.5)/12) im=im+1 TmaxAvgM[im]=TmaxAvgM[im]+Tmax[k] TminAvgM[im]=TminAvgM[im]+Tmin[k] precipAvgM[im]=precip[im]+precip[k] snowAvgM[im]=snow[im]+snow[k] n=n+1 TmaxAvgM[im]=TmaxAvgM[im]/n TminAvgM[im]=TminAvgM[im]/n #Taking an average or sum of a year dateAvgY.append(year[0]+(month[0]-.5)/12) n=0 #day counter iy=0 #year counter for k in range(len(year)): #loop over days if year[k]!=year[k-1] and k>0: #test for new year TmaxAvgY[iy]=TmaxAvgY[iy]/n TminAvgY[iy]=TminAvgY[iy]/n n=0 TmaxAvgY.append(0.) TminAvgY.append(0.) precipAvgY.append(0.) snowAvgY.append(0.) dateAvgY.append(year[k]) if precipAvgY[iy] ==0: precipAvgY[iy]=precipAvgY[iy-1] iy=iy+1 TmaxAvgY[iy]=TmaxAvgY[iy]+Tmax[k] TminAvgY[iy]=TminAvgY[iy]+Tmin[k] precipAvgY[iy]=precipAvgY[iy]+precip[k] snowAvgY[iy]=snowAvgY[iy]+snow[k] n=n+1 TmaxAvgY[iy]=TmaxAvgY[iy]/n TminAvgY[iy]=TminAvgY[iy]/n def C2Flist(c): #function to turn a list from celcius to farenhiet f=[] for k in range(len(c)): f.append(9./5.*c[k]+32.) return f Tmaxf=C2Flist(Tmax) Tminf=C2Flist(Tmin) TmaxAMf=C2Flist(TmaxAvgM) TminAMf=C2Flist(TminAvgM) TmaxAYf=C2Flist(TmaxAvgY) TminAYf=C2Flist(TminAvgY) def AvgList(x): #function to average a list xAvg=0. for k in range(len(x)): xAvg=xAvg+x[k] xAvg=xAvg/(k+1) return xAvg AdateAvgY=array(dateAvgY) ATmaxAYf=array(TmaxAYf) ATminAYf=array(TminAYf) AprecipY=array(precipAvgY) AsnowY=array(snowAvgY) m,b=polyfit(AdateAvgY,ATmaxAYf,1) yfitMax=m*AdateAvgY+b TmaxfPERdec=m*10 TmaxfPERdec=round(TmaxfPERdec,2) print 'Tmax change, deg F/decade', TmaxfPERdec m,b=polyfit(AdateAvgY,ATminAYf,1) yfitMin=m*AdateAvgY+b TminfPERdec=m*10 TminfPERdec=round(TminfPERdec,2) print 'Tmin change, deg F/decade', TminfPERdec #plot(date,Tmaxf,date,Tminf) plot(dateAvgY,TmaxAYf,'r-*',dateAvgY,TminAYf,'b-*', AdateAvgY,yfitMax,'g-',AdateAvgY,yfitMin,'g-') xlabel('Year') ylabel('Temperature Annual Average, f') title(city+': Tmax: '+str(TmaxfPERdec)+' F/decade, Tmin:'+str(TminfPERdec)+' F/decade') legend(['Max Temperature','Min Temperature','Best Fit Line Max', 'Best Fit Line Min'],'center left') # axis([1913,2020,35,70]) grid(True) savefig('figures/'+city+'_Temperature.png') clf() m,b=polyfit(AdateAvgY,AprecipY,1) yfitPrecip=m*AdateAvgY+b precipPERdec=m*100 precipPERdec=round(precipPERdec,2) print 'Precipitation change, deg F/decade', precipPERdec plot(dateAvgY,precipAvgY,'r-*',AdateAvgY,yfitPrecip) xlabel('Year') ylabel('Precipitation Annual Average, inches') title(city+': precip: '+str(precipPERdec)+' inches/century') legend(['Yearly Precipitation','Best Fit Line Precipitation',],'lower left') # axis([1913,2020,35,70]) grid(True) savefig('figures/'+city+'_precip.png') clf() m,b=polyfit(AdateAvgY,AsnowY,1) yfitsnow=m*AdateAvgY+b snowPERdec=m*100 snowPERdec=round(snowPERdec,2) print 'snow change, deg F/decade', snowPERdec plot(dateAvgY,snowAvgY,'r-*',AdateAvgY,yfitsnow) xlabel('Year') ylabel('Snow Annual Sum, inches') title(city+': Snow: '+str(snowPERdec)+' inches/century') legend(['Yearly Snow','Best Fit Line snow',],'upper right') # axis([1913,2020,35,70]) grid(True) savefig('figures/'+city+'_snow.png') clf()