Demo entry 6732859

assignment

   

Submitted by anonymous on Apr 16, 2018 at 02:01
Language: Python. Code size: 10.0 kB.

# -*- coding: utf-8 -*-
"""
Created on Wed Mar 14 21:57:10 2018

@author: TZQ88888
"""
from __future__ import division
import matplotlib.pyplot as plt
import numpy as np
import numpy.ma as ma
import datetime as dt
from scipy.stats.stats import linregress
from netCDF4 import num2date , date2num , Dataset as nc
from mpl_toolkits.basemap import Basemap
import scipy.stats

inputfile_1='MTMG06_ass3_slp.nc'
inputfile_2='MTMG06_ass3_precip.nc'
inputfile_3='MTMG06_ass3_temp.nc'
nc_1=nc(inputfile_1)
#nc_2=nc(inputfile_2)
#nc_3=nc(inputfile_3)


#1.1
msl=nc_1.variables['msl'][:]
lons=nc_1.variables['lon'][:]
lats=nc_1.variables['lat'][:]
raw_times=nc_1.variables['time'][:]
nc_1.close
times=num2date(raw_times,units=nc_1.variables['time'].units,calendar=nc_1.variables['time'].calendar ) 
mean_msl=msl.mean(axis=0)

#plot in basemap
def plot(data,lons=lons, lats=lats,L=10,width=9000000, height =5500000):
    plt.figure()
    
    m= Basemap(width=width,height=height,projection='aeqd',lat_0=50,lon_0=-5,resolution='l')
    m.drawcoastlines()
    X,Y = np.meshgrid(lons, lats) 
    x,y = m(X,Y)
    m.drawmapboundary()
    m.drawparallels(np.arange(0.,90,20.),labels=[1,1,1,1],fontsize=7)
    m.drawmeridians(np.arange(-60.,120.,30.),labels=[1,1,1,1],fontsize=7)

    clevs=np.arange(np.min(data),np.max(data),(np.max(data)-np.min(data))/20)
    cs=m.contourf(x,y,data,clevs,extend='both',cmap='RdYlBu_r' )
    m.colorbar()
    plt.contourf(lons,lats,mean_msl,clevs,extend='both',color='none',hatches=['.' , None]  )
    c=m.contour(x, y,data,L, colors='black', linewidth=.2)
    plt.clabel(c, inline=True, fontsize=10)
    plt.title('Sea level pressure\n',loc='center')
    
plot(mean_msl)
#1.2 PCA
plt.figure()
X=msl
sclats1 = np.sqrt(np.cos(lats*np.pi / 180) )
b=np.swapaxes(X,1,2)
c=b*sclats1
matrix1=np.swapaxes(c,1,2)         
matrix1=np.resize(matrix1,(6882,931))
matrix=matrix1.T
average=np.mean(matrix1,axis=0)
m, n = np.shape(matrix1)
data_adjust = []
avgs = np.tile(average, (m, 1))
data_adjust = matrix1 - avgs

cov = np.cov(data_adjust.T) 
lam,E=np.linalg.eig((cov))
sortlam=np.sort(lam)   #sorf lam 
sortlam=sortlam[-1::-1]
index = np.argsort(-lam)

#PC loading
for i in range(4):
    pc4= np.matrix(E.T[index[i]])*np.sqrt(sortlam[i])

    pc_loading=np.resize(pc4,(19,49))
    b=np.swapaxes(pc_loading,0,1)
    c=b/sclats1
    data_pc4_recon=np.swapaxes(c,0,1)
    plot(data_pc4_recon)         

#1.4
snao= np.matrix(E.T[index[:1]])
snao_data=np.dot(data_adjust,snao.T)
reconData = (snao_data * snao) + average  
plt.figure()
plt.plot(times,snao_data)
plt.xlabel('year')
plt.ylabel('SNAO score')
slope=linregress(np.arange(len(snao_data)),snao_data.T)[0]            
  
#1.4.2
day=np.zeros(shape=(1,19*49))

for i in [0,(np.sqrt(sortlam[0])),(-np.sqrt(sortlam[0]))]:
    day[0,0]=i
    reconData = (day.T * E) + average  
    DATA=np.resize(reconData,(19,49))
    b=np.swapaxes(DATA,0,1)
    c=b/sclats1
    DATA=np.swapaxes(c,0,1)
    plot(DATA,L=1)     


#1.5
lam12= sortlam[:12]
tot =  sum(sortlam)
var_lam= [(i / tot) for i in lam12]
std_lam = [np.sqrt(i)  for i in lam12]              
cum_var_exp = np.cumsum(var_lam)



fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.bar(np.arange(1,13),std_lam,width=0.2)

ax1.set_ylabel('stdev=sqrt(lamda),degrees')
ax1.set_xlabel('PC number')
plt.ylim(0,80)
ax2 = ax1.twinx()
ax2.plot(np.arange(1,13),cum_var_exp,drawstyle= 'steps',color='red')
ax2.set_ylabel('cumulative_variance_explained')
ax2.set_title("first 12 PC")
import scipy.stats as sps
N=6882
lam12=lam12
a=np.sqrt(N/2)
interval=np.arange(24,dtype='float64').reshape(12, 2)
val=sps.norm.interval(0.95,scale=1/np.sqrt(N-1))
for i in range(12):
    interval[i,0]=(val[0]/a*lam12[i]+lam12[i])
    interval[i,1]=(val[1]/a*lam12[i]+lam12[i])

for i in range(12):

    ax1.text(i+0.9,std_lam[i]+2.5, r"%.0f"% (var_lam[i]*100)+'%')
    std_lam_up = [np.sqrt(interval[i,1]) ]
    std_lam_down = [np.sqrt(interval[i,0]) ] 
    ax1.plot([i+0.75,i+1.25],[std_lam_up,std_lam_up],color='black',markersize=6)
    ax1.plot([i+0.75,i+1.25],[std_lam_down,std_lam_down],color='black',markersize=6)

plt.show()

#1.6
d_flood=dt.datetime(2007,7,20,12)
fli= np. where( times == d_flood ) [ 0 ] [ 0 ] 
data=msl[fli]
b=np.swapaxes(data,0,1)
c=b*sclats1
data1=np.swapaxes(c,0,1)
mean=np.mean(data1)
data1=np.resize(data1,19*49)-mean  
for i in [5,20,200]:
    lam=np.matrix(E.T[index[:i]])
    data2=np.dot(data1.T,lam.T)
    reconData = (data2 * lam) + mean
    data3=np.resize(reconData,(19,49))
    b=np.swapaxes(data3,0,1)
    c=b/sclats1
    data4=np.swapaxes(c,0,1)
    plot(data4,width=5000000, height =3500000,L=3)
    



nc_2=nc(inputfile_2)

pre=nc_2.variables['precip'][:]
lons_pre=nc_2.variables['lon'][:]
lats_pre=nc_2.variables['lat'][:]
raw_times_pre=nc_2.variables['time'][:]
nc_2.close
times_pre=num2date(raw_times_pre,units=nc_2.variables['time'].units,calendar=nc_2.variables['time'].calendar )
mean_pre=pre.mean(axis=0) 
def plot_pre(data,lons=lons_pre, lats=lats_pre):
    plt.figure()
    
    m= Basemap(width=5000000, height =3500000,projection='aeqd',lat_0=50,lon_0=15,resolution='l')
    m.drawcoastlines()
    X,Y = np.meshgrid(lons, lats) 
    x,y = m(X,Y)
    m.drawmapboundary()
    m.drawparallels(np.arange(0.,90,20.),labels=[1,1,1,1],fontsize=7)
    m.drawmeridians(np.arange(-60.,120.,30.),labels=[1,1,1,1],fontsize=7)

    clevs=np.arange(np.min(data),np.max(data),0.1)
    cs=m.contourf(x,y,data,clevs,extend='both',cmap='rainbow' )
    m.colorbar()
    #plt.contourf(lons,lats,mean_msl,clevs,extend='both',color='none',hatches=['.' , None]  )
    c=m.contour(x, y,data,10, colors='black', linewidth=.5)
    plt.clabel(c, inline=True, fontsize=10)
    
plot_pre(mean_pre)
A=np.where(snao_data>0)[0] 
B=np.where(snao_data<0)[0] 
C=([np. where( times_pre == i ) [ 0 ] [ 0 ] for i in times[A] if i in times_pre ])
D=([np. where( times_pre == i ) [ 0 ] [ 0 ] for i in times[B] if i in times_pre ])
pos_data=pre[C]   
neg_data=pre[D]
pos_mean=np.mean(pos_data)
neg_mean=np.mean(neg_data)
plot_pre(pos_data.mean(axis=0))
plot_pre(neg_data.mean(axis=0))




def plot_pre1(data,interval,lons=lons_pre, lats=lats_pre):
    plt.figure()
    
    m= Basemap(width=5000000, height =3500000,projection='aeqd',lat_0=50,lon_0=15,resolution='l')
    m.drawcoastlines()
    X,Y = np.meshgrid(lons, lats) 
    x,y = m(X,Y)
    m.drawmapboundary()
    m.drawparallels(np.arange(0.,90,20.),labels=[1,1,1,1],fontsize=7)
    m.drawmeridians(np.arange(-60.,120.,30.),labels=[1,1,1,1],fontsize=7)
    pv_levs=np.arange(-interval,interval,0.05)
    cs=m.contourf(x,y,data,pv_levs,extend='neither',cmap='rainbow'  )
    m.colorbar()
    #plt.contourf(lons,lats,mean_msl,clevs,extend='both',color='none',hatches=['.' , None]  )
    c=m.contour(x, y,data,10, colors='black', linewidth=.5)
    plt.clabel(c, inline=True, fontsize=10)
    

mean=np.mean(pos_data)
std=np.std(pos_data)
interval=sps.norm.interval(0.95,loc=0,scale=std/np.sqrt(N-1))

data_pre=pos_data.mean(axis=0)-neg_data.mean(axis=0)
plot_pre1(data_pre,interval[1])


nc_3=nc(inputfile_3)

temp=nc_3.variables['temp'][:]
lons_temp=nc_3.variables['lon'][:]
lats_temp=nc_3.variables['lat'][:]
raw_times_temp=nc_2.variables['time'][:]
nc_3.close
times_temp=num2date(raw_times_temp,units=nc_3.variables['time'].units,calendar=nc_3.variables['time'].calendar )
mean_temp=temp.mean(axis=0) 

def plot_temp(data,lons=lons_pre, lats=lats_pre):
    plt.figure()
    
    m= Basemap(width=5000000, height =3500000,projection='aeqd',lat_0=50,lon_0=15,resolution='l')
    m.drawcoastlines()
    X,Y = np.meshgrid(lons, lats) 
    x,y = m(X,Y)
    m.drawmapboundary()
    m.drawparallels(np.arange(0.,90,20.),labels=[1,1,1,1],fontsize=7)
    m.drawmeridians(np.arange(-60.,120.,30.),labels=[1,1,1,1],fontsize=7)

    clevs=np.arange(np.min(data),np.max(data),0.1)
    cs=m.contourf(x,y,data,clevs,extend='both',cmap='rainbow' )
    m.colorbar()
    #plt.contourf(lons,lats,mean_msl,clevs,extend='both',color='none',hatches=['.' , None]  )
    c=m.contour(x, y,data,10, colors='black', linewidth=.5)
    plt.clabel(c, inline=True, fontsize=10)
plot_temp(mean_temp)

C=([np. where( times_temp == i ) [ 0 ] [ 0 ] for i in times[A] if i in times_temp ])
D=([np. where( times_temp == i ) [ 0 ] [ 0 ] for i in times[B] if i in times_temp ])
pos_data_temp=temp[C]   
neg_data_temp=temp[D]
pos_mean=np.mean(pos_data_temp)
neg_mean=np.mean(neg_data_temp)
plot_temp(pos_data_temp.mean(axis=0))

N=pos_data_temp.shape[0]
def plot_temp1(data,interval,lons=lons_pre, lats=lats_pre):
    plt.figure()
    
    m= Basemap(width=5000000, height =3500000,projection='aeqd',lat_0=50,lon_0=15,resolution='l')
    m.drawcoastlines()
    X,Y = np.meshgrid(lons, lats) 
    x,y = m(X,Y)
    m.drawmapboundary()
    m.drawparallels(np.arange(0.,90,20.),labels=[1,1,1,1],fontsize=7)
    m.drawmeridians(np.arange(-60.,120.,30.),labels=[1,1,1,1],fontsize=7)
    pv_levs=np.arange(-interval,interval,0.01)
    cs=m.contourf(x,y,data,pv_levs,extend='neither',cmap='rainbow'  )
    m.colorbar()
    #plt.contourf(lons,lats,mean_msl,clevs,extend='both',color='none',hatches=['.' , None]  )
    c=m.contour(x, y,data,10, colors='black', linewidth=.5)
    plt.clabel(c, inline=True, fontsize=10)
    
mean=np.mean(pos_data_temp)
std=np.std(pos_data_temp)
interval=sps.norm.interval(0.95,loc=0,scale=std/np.sqrt(N-1))

data_temp=pos_data_temp.mean(axis=0)-neg_data_temp.mean(axis=0)
plot_temp1(data_temp,interval[1])



tol_temp_n=0
temp_n=0

for i in range(36):
    for k in range(72):
        if abs(data_temp[i][k])<interval[1]:
            temp_n+=1
        if abs(data_temp[i][k])<1000:
            tol_temp_n+=1


dif_grid_num=tol_temp_n-temp_n

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).