1 R語言
1.1 泰勒圖
1.1.1 plotrix包
library(plotrix)ref<-rnorm(30,sd=2)model1<-ref+rnorm(30)/2model2<-ref+rnorm(30)oldpar<-taylor.diagram(ref,model1)taylor.diagram(ref,model2,add=TRUE,col="blue")lpos<-1.5*sd(ref)legend(lpos,lpos,legend=c("Better","Worse"),pch=19,col=c("red","blue"))par(oldpar)taylor.diagram(ref,model1,pos.cor=FALSE)taylor.diagram(ref,model2,add=TRUE,col="blue")
1.1.2 openair包
library(openair)dat <- selectByDate(mydata, year = 2003)dat <- data.frame(date = mydata$date, obs = mydata$nox, mod = mydata$nox)## now make mod worse by adding bias and noise according to the month## do this for 3 different modelsdat <- transform(dat, month = as.numeric(format(date, "%m")))mod1 <- transform(dat, mod = mod + 10 * month + 10 * month * rnorm(nrow(dat)),model = "model 1")## lag the results for mod1 to make the correlation coefficient worse## without affecting the sdmod1 <- transform(mod1, mod = c(mod[5:length(mod)], mod[(length(mod) - 3) :length(mod)]))## model 2mod2 <- transform(dat, mod = mod + 7 * month + 7 * month * rnorm(nrow(dat)),model = "model 2")## model 3mod3 <- transform(dat, mod = mod + 3 * month + 3 * month * rnorm(nrow(dat)),model = "model 3")mod.dat <- rbind(mod1, mod2, mod3)## basic Taylor plotTaylorDiagram(mod.dat, obs = "obs", mod = "mod", group = "model")## Taylor plot by seasonTaylorDiagram(mod.dat, obs = "obs", mod = "mod", group = "model", type = "season")## now show how to evaluate model improvement (or otherwise)mod1a <- transform(dat, mod = mod + 2 * month + 2 * month * rnorm(nrow(dat)),model = "model 1")mod2a <- transform(mod2, mod = mod * 1.3)mod3a <- transform(dat, mod = mod + 10 * month + 10 * month * rnorm(nrow(dat)),model = "model 3")mod.dat2 <- rbind(mod1a, mod2a, mod3a)mod.dat$mod2 <- mod.dat2$mod## now we have a data frame with 3 models, 1 set of observations## and TWO sets of model predictions (mod and mod2)## do for all modelsTaylorDiagram(mod.dat, obs = "obs", mod = c("mod", "mod2"), group = "model")## End(Not run)## Not run:## all models, by seasonTaylorDiagram(mod.dat, obs = "obs", mod = c("mod", "mod2"), group = "model",type = "season")## consider two groups (model/month). In this case all months are shown by model## but are only differentiated by model.TaylorDiagram(mod.dat, obs = "obs", mod = "mod", group = c("model", "month"))1.2 bcp包--貝葉斯突變檢驗
library(bcp)##### univariate sequential data ###### an easy problem with 2 true change pointsset.seed(5)x <- c(rnorm(50), rnorm(50, 5, 1), rnorm(50))bcp.1a <- bcp(x)plot(bcp.1a, main="Univariate Change Point Example")legacyplot(bcp.1a)# a hard problem with 1 true change pointset.seed(5)x <- rep(c(0,1), each=50)y <- x + rnorm(50, sd=1)bcp.1b <- bcp(y)plot(bcp.1b, main="Univariate Change Point Example")##### multivariate sequential data ###### an easy problem in k=3 dimensionsset.seed(5)x <- rnorm(6, sd=3)y <- rbind(cbind(rnorm(50, x[1]), rnorm(50, x[2]), rnorm(50, x[3])),cbind(rnorm(50, x[4]), rnorm(50, x[5]), rnorm(50, x[6])))bcp.2a <- bcp(y)plot(bcp.2a, main="Multivariate (k=3) Change Point Example")plot(bcp.2a, separated=TRUE, main="Multivariate (k=3) Change Point Example")# a harder problem in k=5 dimensionsset.seed(5)means1 <- rep(0, 5)means2 <- rep(1, 5)x <- rbind(matrix(rep(means1, each=50), nrow=50),matrix(rep(means2, each=50), nrow=50))y <- x + rnorm(length(x), sd=1)bcp.2b <- bcp(cbind(y))plot(bcp.2b, main="Multivariate (k=5) Change Point Example")##### linear models with sequential data ###### 1 true change point at location 50; the predicting variable x is not related to locationx <- rnorm(100)b <- rep(c(3,-3), each=50)y <- b*x + rnorm(100)bcp.3a <- bcp(y, x)# in the two plots that follow, the location IDs are used as the plot characterspar(mfrow=c(1,2))plot(y ~ x, type="n", main="Linear Regression: Raw Data")text(x, y, as.character(1:100), col=(b/3)+2)plot(y ~ x, type="n", main="Linear Regression: Posterior Means")text(x, bcp.3a$posterior.mean[,1], as.character(1:100), col=(b/3)+2)plot(bcp.3a, main="Linear Regression Change Point Example")# 1 true change point at location 50; the predicting variable x is equal to locationx <- 1:100b <- rep(c(3,-3), each=50)y <- b*x + rnorm(100, sd=50)bcp.3b <- bcp(y, x)plot(bcp.3b, main="Linear Regression Change Point Example")##### univariate data on a grid ####### Not run:set.seed(5)adj <- makeAdjGrid(20)z <- rep(c(0, 2), each=200)y <- z + rnorm(400, sd=1)out <- bcp(y, adj=adj, burnin=500, mcmc=500)if (require("ggplot2")) {df <- data.frame(mean=z, data = y, post.means = out$posterior.mean[,1],post.probs = out$posterior.prob,i = rep(1:20, each=20), j = rep(1:20, times=20))# visualize the meansg <- ggplot(df, aes(i,j)) +geom_tile(aes(fill = mean), color='white') +scale_fill_gradientn(limits=range(y), colours=c('white', 'steelblue'))+ggtitle("True Means")print(g)# visualize the datag <- ggplot(df, aes(i,j)) +geom_tile(aes(fill = data), color='white') +scale_fill_gradientn(limits=range(y), colours=c('white', 'steelblue'))+ggtitle("Observed Data")print(g)# visualize the posterior means/probsg <- ggplot(df, aes(i,j)) +geom_tile(aes(fill = post.means), color='white') +scale_fill_gradientn(limits=range(y), colours=c('white', 'steelblue'))+ggtitle("Posterior Means")print(g)g <- ggplot(df, aes(i,j)) +geom_tile(aes(fill = post.probs), color='white') +scale_fill_gradientn(limits=c(0, 1), colours=c('white', 'steelblue'))+ggtitle("Posterior Boundary Probabilities")print(g)}## End(Not run)## Not run:##### multivariate data on a grid #####set.seed(5)x <- rnorm(6, sd=3)y <- rbind(cbind(rnorm(50, x[1]), rnorm(50, x[2]), rnorm(50, x[3])),cbind(rnorm(50, x[4]), rnorm(50, x[5]), rnorm(50, x[6])))adj <- makeAdjGrid(10)a <- bcp(y, adj=adj, p0=0.4, burnin=500, mcmc=500)##### linear models on a grid #####set.seed(5)x <- rnorm(100)b <- rep(c(3,-3), each=50)y <- b*x + rnorm(100)adj <- makeAdjGrid(10)a <- bcp(y,x,adj=adj, p0=0.4, burnin=500, mcmc=500)##### linear models on a grid using pseudo-APPs #####x <- rnorm(100)b <- rep(c(3,-3), each=50)y <- b*x + rnorm(100)adj <- makeAdjGrid(10)a <- bcp(y,x,adj=adj, p0=0.4, burnin=500, mcmc=500, p1 = 0)## End(Not run)##### univariate data on a graph ####### Not run:demo(bcpgraph)## End(Not run)###### Real Data Examples ######## Not run:# Coriell chromosome 11: univariate sequential datademo(coriell)# U.S. ex-post interest rate: univariate sequential datademo(RealInt)# Lombard: univariate sequential data (with and without linear models)demo(Lombard)# Quebec rivers: multivariate sequential datademo(QuebecRivers)# New Haven housing: linear models on a graphdemo(NewHaven)## End(Not run)1.3 visreg包
library(visreg)par(mfrow = c(2, 2)) ##reg01模型中有四個變量，設置為2 x 2的圖visreg(reg01)visreg(reg01, "x1") ##注意，變量名一定要用引號括起來，不然會報錯visreg(reg01, "cat")visreg(reg01, "x1", by = "cat")visreg(reg01, "x1", by = "x2", breaks = 4, layout = c(2, 2))visreg(reg01, "x1", by = "cat", overlay = TRUE)
1.4 ggcor包
devtools::install_github("houyunhuang/ggcor") library(devtools)devtools::install_github("houyunhuang/ggcor", force = TRUE) library(ggcor)
2 Matlab
2.1 GreenSim團隊
%時間序列的貝葉斯突變檢測算法MATLAB源碼function P=BMDCP(Y,R,M)%% Bayesian Method for Detecting Change Points%% 輸入參數列表% Y 時間序列原始數據% R 控制參數Var0/Var，Var為全序列方差，Var0為均值方差，R應大於等於4% M 端點間隔，取值3～5%% 輸入參數列表% P 突變點的概率分布%% 第一步：數據歸一化處理%歸一化處理的目的在於避免浮點計算時溢出，屬程序設計時考慮的問題，它不影響結果Y=Y./max(Y);%% 第二步：計算後驗分布的數字特徵N=length(Y);%序列的長度Mu0=mean(Y);%全序列的均值Var=var(Y);%全序列方差MuA=zeros(N,1);MuB=zeros(N,1);VarB=zeros(N,1);for i=M:(N-M) YA=Y(1:i);%左子序列 YB=Y((i+1):N);%右子序列 VarA(i)=Var/(1/R+i); VarB(i)=Var/(1/R+N-i);endMuA(1:(M-1))=MuA(M);MuB(1:(M-1))=MuB(M);MuA((N-M+1):N)=MuA(N-M);MuB((N-M+1):N)=MuB(N-M);%% 第三步：計算後驗概率分布P=zeros(N,1);for k=M:(N-M) A=0; for i=1:k % 注意：這裡對似然函數進行了取對數處理，也是為了防止CPU溢出 a=log((1/sqrt(2*pi*Var))*exp(-(Y(i)-MuA(i))^2/(2*Var))); end B=0; for j=(k+1):N b=log((1/sqrt(2*pi*Var))*exp(-(Y(j)-MuB(j))^2/(2*Var))); end P(k)=A+B;end%% 第四步：計算突變點的概率分布圖P=exp(P);P=P/sum(P);maxP=max(P);bar(1:N,P);axis([1,N,0,1.2*maxP]);hold onp0=(1/N)*ones(N,1);plot(1:N,p0,'--r');legend('後驗概率','先驗概率');xlabel('k (month)','FontName','TimesNewRoman','FontSize',12);ylabel('P(k)','FontName','TimesNewRoman','Fontsize',12);%非平穩時間序列突變檢測的啟發式分割算法（BG算法）function [FLAG,AllT,AllTmax,AllPTmax]=BGA(X,P0,L0)%% 非平穩時間序列突變檢測的啟發式分割算法——BG算法
%% 輸入參數列表% X 待檢測的數據，列向量存儲% P0 顯著性水平門限值，低於此值的不再分割% L0 最小分割尺度，子段長度小於此值的不再分割%% 輸出參數列表% FLAG 分割點標記，列向量存儲，長度與X相同% AllT 與分割點對應的全部t檢驗序列，其首位數字為起點坐標% AllTmax 與分割點對應的全部t檢驗序列的最大值% AllPTmax 與分割點對應的全部t檢驗序列對應的統計顯著性%% 第一步：變量初始化N=length(X);FLAG=zeros(N,1);FLAG(1)=0.1;FLAG(N)=0.1;AllT=cell(0,0);AllTmax=cell(0,0);AllPTmax=cell(0,0);%% 第二步：產生第一個突變點，並對序列進行分割[T,Tmax,p,PTmax]=Tseries(X);T=[1;T];if PTmax<P0 flag=FLAG; flag(1)=0; flag(N)=0; pos3=flag(pos2); returnend%記錄輸出數據FLAG(p)=1;AllT=[AllT;T];AllTmax=[AllTmax;Tmax];AllPTmax=[AllPTmax;PTmax];%以下為兩個控制計數器counter=2;%下一個突變點的序號TC=0;%臨時計數器%%while 1%設置死循環%% 第三步：對每個段進行突變檢測，能分割則分割，直到不能分割為止 pos=find(FLAG>0); M=length(pos)-1;%當前子段數目 for m=1:M s=pos(m); t=pos(m+1); L=length(SubX); if L>=L0 [T,Tmax,p,PTmax]=Tseries(SubX); T=[s;T]; if PTmax>=P0 TC=TC+1; FLAG(s+p-1)=counter; AllT=[AllT;T]; AllPTmax=[AllPTmax;PTmax]; counter=counter+1; end end end%% 第四步：返回輸出數據 if TC==0 flag=FLAG; flag(1)=0; flag(N)=0; pos3=flag(pos2); FLAG=[pos2,pos3]; return end%% TC=0;%%endfunction [T,Tmax,p,PTmax]=Tseries(x)%% 計算t檢驗統計序列的子函數
%% 參數列表% x 時間序列，N×1列向量% T t檢驗序列，N×1列向量% Tmax t檢驗序列的最大值% p t檢驗序列最大值對應的下標% PTmax Tmax對應的統計顯著性%% 參數初始化N=length(x);T=zeros(N,1);%% 以下是主循環，用於創建t檢驗序列for i=3:(N-2)%最左邊以及最右邊的兩個點沒有對應的t檢驗值（或者說，其值初始化為0） x1=x(1:i);%序列左邊部分 N1=length(x1);%左邊序列的長度 x2=x(i:N);%序列右邊部分 N2=length(x2);%右邊序列的長度 mean_x2=mean(x2);%右邊部分的均值 std_x2=std(x2);%右邊部分的標準差 %下面是計算合併偏差的公式，中英文文獻裡的這個公式略有不同，此處以英文文獻為準 SD=sqrt(1/N1+1/N2)*sqrt(((N1-1)*std_x1^2+(N2-1)*std_x2^2)/(N1+N2-2)); T(i)=abs((mean_x1-mean_x2)/SD);end%% 計算其它三個輸出參數Tmax=max(T);%t檢驗序列的最大值pos=find(T==Tmax);Eta=4.19*log(N)-11.54;%計算PTmax用的參數Delta=0.40;%計算PTmax用的參數v=N-2;%計算PTmax用的參數c=v/(v+Tmax^2);%不完全B函數的下標PTmax=(1-betainc(c,Delta*Eta,Delta));%調用不完全beta函數%MannKendall突變檢測算法MATLAB源碼%% 以下是單邊MannKendall趨勢檢測算法clc;clear;close allload dataX=MJL;N=length(X);U=zeros(N-1,1);for t=2:N x=X(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n S=S+sign(x(j)-x(k)); end end VarS=n*(n-1)*(2*n+5)/18; if S>0 Z=(S+1)/sqrt(VarS); elseif S==0 Z=0; else Z=(S-1)/sqrt(VarS); end U(t-1)=Z;endfigure(1)plot(1:(N-1),U,'linewidth',1.5);hold onplot(1:(N-1),1.96*ones(N-1,1),':','linewidth',1);legend('統計量','0.05顯著水平');hold onplot(1:(N-1),0*ones(N-1,1),'-.','linewidth',1);plot(1:(N-1),1.96*ones(N-1,1),':','linewidth',1);plot(1:(N-1),-1.96*ones(N-1,1),':','linewidth',1);axis([1,N-1,-5,5]);xlabel('t (month)','FontName','TimesNewRoman','FontSize',12);ylabel('U統計量','FontName','TimesNewRoman','Fontsize',12);%grid on%計算趨勢顯著水平Alpha=1-normcdf(U(end),0,1);disp('趨勢的顯著水平為');disp(Alpha);%計算整體趨勢的變化速率Qi=zeros(N*(N-1)/2,1);counter=1;for k=1:(N-1) for j=(k+1):N Qi(counter)=(X(j)-X(k))/(j-k); counter=counter+1; endendQ=median(Qi);disp('整體趨勢的變化速率為');disp(Q);
%% 以下是雙邊MannKendall突變檢測算法clc;clear;close allload dataX=YJL;%計算UF統計量N=length(X);UF=zeros(N-1,1);for t=2:N x=X(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n if x(j)>x(k) S=S+1; else S=S+0; end end end ES=n*(n+1)/4; VarS=n*(n-1)*(2*n+5)/72; Z=(S-ES)/sqrt(VarS); UF(t-1)=Z;end%計算UB統計量Y=flipud(X);UB=zeros(N-1,1);for t=2:N x=Y(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n if x(j)>x(k) S=S+1; else S=S+0; end end end ES=n*(n+1)/4; VarS=n*(n-1)*(2*n+5)/72; Z=(S-ES)/sqrt(VarS); UB(t-1)=-Z;end%繪圖figure(2)plot(1:(N-1),UF,'r-','linewidth',1.5);hold onplot(1:(N-1),UB,'b-.','linewidth',1.5);plot(1:(N-1),1.96*ones(N-1,1),':','linewidth',1);axis([1,N-1,-4,4]);legend('UF統計量','UB統計量','0.05顯著水平');xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);ylabel('統計量','FontName','TimesNewRoman','Fontsize',12);%grid onhold onplot(1:(N-1),0*ones(N-1,1),'-.','linewidth',1);plot(1:(N-1),1.96*ones(N-1,1),':','linewidth',1);plot(1:(N-1),-1.96*ones(N-1,1),':','linewidth',1);%水電站中長期優化調度的粒子群算法MATLAB源碼%% 參數設置A=8.5;%出力係數，常數Tt=730*ones(12,1);%第t個時段的小時數%注意：一年按363天*24小時算，均分為12個月HtLB=55*ones(12,1);%第t時段水位約束的下界，單位：米HtUB=[65;65;65;61;61;61;61;65;65;65;65;65];%第t時段水位約束的上界，單位：米VtLB=zeros(12,1);VtUB=zeros(12,1);for i=1:12 VtLB(i)=Ht2Vt(HtLB(i)); VtUB(i)=Ht2Vt(HtUB(i));end%注意：蓄水量Vt和水位Ht之間有一一對應的關係，單位：立方米NtLB=260000*ones(12,1);%出力約束的下界，單位：千瓦NtUB=1400000*ones(12,1);%出力約束的上屆，單位：千瓦%注意：Nt=A*Qt*HtQtLB=308*ones(12,1);%洩流量下界，單位：立方米/秒QtUB=29200*ones(12,1);%洩流量上界，單位：立方米/秒qt=[373;859;1568;2100;3210;5049;1596;1160;925;781;572;1010];%入庫流量，單位：立方米/秒%注意：以上三個量，時間單位相乘時，小時乘以3600轉化成秒%% 調用粒子群算法K=60;N=80;w=0.5;c1=0.3;c2=0.2;[BESTX,BESTY,ALLX,ALLY]=PSO(K,N,w,c1,c2,VtLB,VtUB,QtLB,QtUB,NtLB,NtUB,qt,A,Tt);%%X=BESTX{K};[Vt,Qt,St]=DeCode(X);disp('最佳蓄水量');disp(Vt);disp('最佳平均發電流量');disp(Qt);disp('最佳平均棄水流量');disp(St);disp('最大總發電量');disp(-BESTY(K));
function [BESTX,BESTY,ALLX,ALLY]=PSO(K,N,w,c1,c2,VtLB,VtUB,QtLB,QtUB,NtLB,NtUB,qt,A,Tt)%% 此函數實現粒子群優化算法，用於水電站調度優化%% 輸入參數列表% K 迭代次數% N 粒子個數% w 慣性因子% c1 加速因子，針對歷史最優狀態% c2 加速因子，針對全局最優狀態% VtLB 蓄水量約束下界，立方米，N*1向量% VtUB 蓄水量約束上界，立方米，N*1向量% QtLB 平均發電流量約束下界，立方米/秒，N*1向量% QtUB 平均發電流量約束上界，立方米/秒，N*1向量% NtLB 出力約束下界，千瓦，N*1向量% NtUB 出力約束上界，千瓦，N*1向量% qt 入庫流量，立方米/秒，N*1向量% A 出力係數% Tt 一個周期的小時數%% 輸出參數列表% BESTX K×1細胞結構，每一個元素是M×1向量，記錄每一代的最優個體% BESTY K×1矩陣，記錄每一代的最優個體的評價函數值% ALLX K×1細胞結構，每一個元素是M×N矩陣，記錄全部個體% ALLY K×N矩陣，記錄全部個體的評價函數值%% 第一步：粒子狀態初始化M=3*length(VtLB);farm=zeros(N,M);for i=N [Vt,Qt,St]=Initialize2(VtLB,VtUB,QtLB,QtUB,NtLB,NtUB,qt,A,Tt); X=EnCode(Vt,Qt,St); farm(i,:)=X';end%粒子歷史最優狀態初始化Pfarm=farm;%粒子群全局最優狀態for i=1:N X=Pfarm(i,:); [Vt,Qt,St]=DeCode(X'); SE=ObjFun(A,Qt,St,Vt,Tt); Fitness(i)=SE;endMinFit=min(Fitness);POS=find(Fitness==MinFit);Gfarm=Pfarm(POS(1),:);%輸出變量初始化ALLX=cell(K,1);%細胞結構，每一個元素是M×N矩陣，記錄每一代的個體ALLY=zeros(K,N);%K×N矩陣，記錄每一代評價函數值BESTX=cell(K,1);%細胞結構，每一個元素是M×1向量，記錄每一代的最優個體BESTY=zeros(K,1);%K×1矩陣，記錄每一代的最優個體的評價函數值k=1;%迭代計數器初始化%% 第二步：迭代過程while k<=K%% 粒子狀態更新 newfarm=farm; for i=1:N newfarm(i,:)=w*farm(i,:)+c1*rand*(Pfarm(i,:)-farm(i,:))+c2*rand*(Gfarm-farm(i,:)); X=newfarm(i,:); [Vt,Qt,St]=DeCode(X'); [Flag,Vt,Qt,St]=Correct(Vt,Qt,St,VtLB,VtUB,QtLB,QtUB,NtLB,NtUB,qt,A,Tt); X=EnCode(Vt,Qt,St); newfarm(i,:)=X'; end%% 歷史最優狀態和全局最優狀態更新 NEWFIT=zeros(1,N); for i=1:(N) aa=newfarm(i,:); [Vt,Qt,St]=DeCode(aa'); FitA=ObjFun(A,Qt,St,Vt,Tt); bb=Pfarm(i,:); [Vt,Qt,St]=DeCode(bb'); FitB=ObjFun(A,Qt,St,Vt,Tt); NEWFIT(i)=FitA; if FitA<FitB Pfarm(i,:)=A; end if FitA<Gfarm Gfarm=A; end end%% 記錄 ALLX{k}=newfarm; ALLY(k,:)=NEWFIT; BESTX{k}=Gfarm; [Vt,Qt,St]=DeCode(Gfarm'); BESTY(k)=ObjFun(A,Qt,St,Vt,Tt); disp(k); k=k+1;end%基於遺傳算法的投影尋蹤模型Matlab源碼%% 「投影尋蹤＋遺傳算法優化」的主仿真程序%% 第一步：仿真參數設置clearclcclose allload Q5.txtDD=Q5;%導入D矩陣[n,p]=size(DD);np=15; %訓練樣本的個數，前面1～np個樣本用於建立模型，剩下的樣本用於預測if np>=n error('用於預測的樣本個數不能大於或等於樣本總數，請重新設置');endyear=1:np;%選擇參與計算的樣本，默認選擇全部Factor=1:p;%選擇部分指標，默認選擇全部D=DD(year,Factor);K=50; %迭代次數N=30; %種群規模Pm=0.3; %變異概率LB=-ones(1,p); %決策變量的下界UB=ones(1,p); %決策變量的上界Alpha=0.1; %窗口半徑係數，典型取值0.1b%% 調用遺傳算法優化投影尋蹤模型的程序[BESTX,BESTY,ALLX,ALLY]=GAUCP(K,N,Pm,LB,UB,D,Alpha)%% 以下均為整理輸出結果%所有數據都在workspace裡，最值得關注的三個數據是% Z 投影指標值，和參考文獻裡的符號是一致的% Best_a 最佳投影向量，參考文獻裡也是用的符號a，這裡加了個前綴Best，表示最佳% BESTY 投影尋蹤模型中的目標函數的變化情況，文獻中的模型是最大化模型，這裡按照慣例，對其加了個負號成為最小化模型Best_a=(BESTX{K})';%方向向量disp('最佳投影向量為');disp(Best_a);d=zeros(np,p);DDjmax=max(DD);DDjmin=min(DD);for i=1:np d(i,:)=(DD(i,:)-DDjmin)./(DDjmax-DDjmin);endZ=zeros(np,1);for i=1:np Z(i)=abs(sum(Best_a.*d(i,:)));endZ=abs(Z);%%figure(2)%投影散布圖plot(year,abs(Z),'bd','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','b','MarkerSize',5);%axis([1,12,0,2.5]);%圖形邊界根據需要顯示grid onxlabel('Year','FontName','Times New Roman','FontSize',12);ylabel('Projective Value','FontName','Times New Roman','Fontsize',12);%%figure(3)[newZ,I]=sort(Z);newyear=year(I);plot(year,abs(newZ),'bd','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','b','MarkerSize',5);%axis([1,12,0,2.5]);%圖形邊界根據需要顯示grid onxlabel('Year','FontName','Times New Roman','FontSize',12);ylabel('Projective Value','FontName','Times New Roman','Fontsize',12);%%n2=n-np;d2=zeros(n2,p);for i=1:n2 d2(i,:)=(DD(i+np,:)-DDjmin)./(DDjmax-DDjmin);endZ2=zeros(n2,1);for i=1:n2 Z2(i)=abs(sum(Best_a.*d2(i,:)));endZ2=abs(Z2);disp('預測樣本的投影預測值為');disp(Z2);%%figure(4)%投影散布圖plot([Z;Z2],'bd','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','b','MarkerSize',5);hold onplot((np+1):n,Z2,'bo','LineWidth',1,'MarkerEdgeColor','r','MarkerFaceColor','r','MarkerSize',5);legend('訓練樣本投影值','預測樣本投影值');%axis([1,12,0,2.5]);%圖形邊界根據需要顯示grid onxlabel('Year','FontName','Times New Roman','FontSize',12);ylabel('Projective Value','FontName','Times New Roman','Fontsize',12);%支持向量機和BP神經網絡比較研究function [Alpha1,Alpha2,Alpha,Flag,B]=SVMNR(X,Y,Epsilon,C,TKF)%%% SVMNR.m% Support Vector Machine for Nonlinear Regression% All rights reserved%%% 支持向量機非線性回歸通用程序% 程序功能：% 使用支持向量機進行非線性回歸，得到非線性函數y=f(x1,x2,…,xn)的支持向量解析式，% 求解二次規劃時調用了優化工具箱的quadprog函數。本函數在程序入口處對數據進行了% [-1,1]的歸一化處理，所以計算得到的回歸解析式的係數是針對歸一化數據的，仿真測% 試需使用與本函數配套的Regression函數。% 主要參考文獻:% 朱國強,劉士榮等.支持向量機及其在函數逼近中的應用.華東理工大學學報% 輸入參數列表% X 輸入樣本原始數據，n×l的矩陣，n為變量個數，l為樣本個數% Y 輸出樣本原始數據，1×l的矩陣，l為樣本個數% Epsilon ε不敏感損失函數的參數，Epsilon越大，支持向量越少% C 懲罰係數，C過大或過小，泛化能力變差% TKF Type of Kernel Function 核函數類型% TKF=1 線性核函數，注意：使用線性核函數，將進行支持向量機的線性回歸% TKF=2 多項式核函數% TKF=3 徑向基核函數% TKF=4 指數核函數% TKF=5 Sigmoid核函數% TKF=任意其它值，自定義核函數% 輸出參數列表% Alpha1 α係數% Alpha2 α*係數% Alpha 支持向量的加權係數（α－α*）向量% Flag 1×l標記，0對應非支持向量，1對應邊界支持向量，2對應標準支持向量% B 回歸方程中的常數項%----%%%---數據歸一化處理---nntwarn offX=premnmx(X);Y=premnmx(Y);%%%%%---核函數參數初始化-switch TKF case 1 %線性核函數 K=sum(x.*y) %沒有需要定義的參數 case 2 %多項式核函數 K=(sum(x.*y)+c)^p c=0.1; p=2; case 3 %徑向基核函數 K=exp(-(norm(x-y))^2/(2*sigma^2)) sigma=10; case 4 %指數核函數 K=exp(-norm(x-y)/(2*sigma^2)) sigma=10; case 5 %Sigmoid核函數 K=1/(1+exp(-v*sum(x.*y)+c)) v=0.5; c=0; otherwise %自定義核函數，需由用戶自行在函數內部修改，注意要同時修改好幾處！ %暫時定義為 K=exp(-(sum((x-y).^2)/(2*sigma^2))) sigma=8;end%%%%%---構造K矩陣---l=size(X,2);K=zeros(l,l);%K矩陣初始化for i=1:l for j=1:l x=X(:,i); y=X(:,j); switch TKF%根據核函數的類型，使用相應的核函數構造K矩陣 case 1 K(i,j)=sum(x.*y); case 2 K(i,j)=(sum(x.*y)+c)^p; case 3 K(i,j)=exp(-(norm(x-y))^2/(2*sigma^2)); case 4 K(i,j)=exp(-norm(x-y)/(2*sigma^2)); case 5 K(i,j)=1/(1+exp(-v*sum(x.*y)+c)); otherwise K(i,j)=exp(-(sum((x-y).^2)/(2*sigma^2))); end endend%%%%%--構造二次規劃模型的參數H,Ft,Aeq,Beq,lb,ub----%支持向量機非線性回歸，回歸函數的係數，要通過求解一個二次規劃模型得以確定Beq=0;lb=eps.*ones(2*l,1);ub=C*ones(2*l,1);%%%%%----調用優化工具箱quadprog函數求解二次規劃----OPT=optimset;OPT.LargeScale='off';OPT.Display='off';%%%%%----整理輸出回歸方程的係數Alpha1=(Gamma(1:l,1))';Alpha=Alpha1-Alpha2;Flag=2*ones(1,l);%%%%%--支持向量的分類----Err=0.000000000001;for i=1:l AA=Alpha1(i); BB=Alpha2(i); if (abs(AA-0)<=Err)&&(abs(BB-0)<=Err) Flag(i)=0;%非支持向量 end if (AA>Err)&&(AA<C-Err)&&(abs(BB-0)<=Err) Flag(i)=2;%標準支持向量 end if (abs(AA-0)<=Err)&&(BB>Err)&&(BB<C-Err) Flag(i)=2;%標準支持向量 end if (abs(AA-C)<=Err)&&(abs(BB-0)<=Err) Flag(i)=1;%邊界支持向量 end if (abs(AA-0)<=Err)&&(abs(BB-C)<=Err) Flag(i)=1;%邊界支持向量 endend%%%%%計算回歸方程中的常數項B---B=0;counter=0;for i=1:l AA=Alpha1(i); BB=Alpha2(i); if (AA>Err)&&(AA<C-Err)&&(abs(BB-0)<=Err) %計算支持向量加權值 SUM=0; for j=1:l if Flag(j)>0 switch TKF case 1 SUM=SUM+Alpha(j)*sum(X(:,j).*X(:,i)); case 2 SUM=SUM+Alpha(j)*(sum(X(:,j).*X(:,i))+c)^p; case 3 SUM=SUM+Alpha(j)*exp(-(norm(X(:,j)-X(:,i)))^2/(2*sigma^2)); case 4 SUM=SUM+Alpha(j)*exp(-norm(X(:,j)-X(:,i))/(2*sigma^2)); case 5 SUM=SUM+Alpha(j)*1/(1+exp(-v*sum(X(:,j).*X(:,i))+c)); otherwise SUM=SUM+Alpha(j)*exp(-(sum((X(:,j)-X(:,i)).^2)/(2*sigma^2))); end end end B=B+b; counter=counter+1; end if (abs(AA-0)<=Err)&&(BB>Err)&&(BB<C-Err) SUM=0; for j=1:l if Flag(j)>0 switch TKF case 1 SUM=SUM+Alpha(j)*sum(X(:,j).*X(:,i)); case 2 SUM=SUM+Alpha(j)*(sum(X(:,j).*X(:,i))+c)^p; case 3 SUM=SUM+Alpha(j)*exp(-(norm(X(:,j)-X(:,i)))^2/(2*sigma^2)); case 4 SUM=SUM+Alpha(j)*exp(-norm(X(:,j)-X(:,i))/(2*sigma^2)); case 5 SUM=SUM+Alpha(j)*1/(1+exp(-v*sum(X(:,j).*X(:,i))+c)); otherwise SUM=SUM+Alpha(j)*exp(-(sum((X(:,j)-X(:,i)).^2)/(2*sigma^2))); end end end b=Y(i)-SUM+Epsilon; counter=counter+1; endendif counter==0 B=0;else B=B/counter;end2.2 颱風路徑
clc;clear;close allA=xlsread('testdata.xlsx');%%[a1]=find(A(:,1)==66666);j=1;l=1;B=[];C=[];for k=1:length(a1) if 220>A(a1(k)+1,3)&&A(a1(k)+1,3)>50&&1200>A(a1(k)+1,4)&&A(a1(k)+1,4)>1050 B(j,:)=A(a1(k)+1,:);% 南海生成的颱風 j=j+1; else C(l,:)=A(a1(k)+1,:);% 非局地颱風 l=l+1; endend%% lat=B(:,3)/10;lon=B(:,4)/10;lat1=C(:,3)/10;lon1=C(:,4)/10;figurem_proj('Equidistant Cylindrical','lat',[0 40],'long',[100 180]); %投影方式及繪圖範圍設定m_gshhs_l('patch',[0.98 0.98 0.98],'EdgeColor',[0 0 0]); %線條及色塊顏色設定% m_grid(); %邊框及網格設定%m_grid('linestyle','none','fontsize',18);m_grid('box','fancy','tickdir','in');hold onm_plot(lon,lat,'.k');m_plot(lon1,lat1,'.r');m_plot(repmat(120,1,18),5:22,'color','m','linewidth',1.5)m_plot(repmat(105,1,18),5:22,'color','m','linewidth',1.5)m_plot(105:5:120,repmat(22,1,4),'color','m','linewidth',1.5)m_plot(105:5:120,repmat(5,1,4),'color','m','linewidth',1.5)%% 統計兩類颱風月頻數、年頻數D=B; %B是土颱風，即局地颱風數，統計非局地颱風時注釋此行%D=C; %C是非局地颱風數time_D=num2str(D(:,1));n=68;%年數y=1949:2016;Y=reshape(repmat(y,12,1),n*12,1);mm=1:12;M=reshape(repmat(mm',1,n),n*12,1);m=1;d(:,3)=zeros(n*12,1);d(:,1)=Y;d(:,2)=M;%d為颱風頻數統計j=1;for i=1:n*12 k=0; for j=1:length(D(:,1)) if Y(i)==str2num(time_D(j,1:4))&&M(i)==str2num(time_D(j,5:6)) k=k+1; end d(i,3)=k; endend x=1:12;y=1949:2016; D_fre_mon=nanmean(reshape(d(:,3),12,n),2);%平均月頻數 D_fre_yr=sum(reshape(d(:,3),12,n),1)';%颱風年頻數 D_fre_anom = D_fre_yr- mean(D_fre_yr);% 距平%% 月頻數柱狀圖 figurebar(x,D_fre_mon,0.4,'FaceColor',[0.4 0.4 0.4]);set(gca,'FontSize',12);ylabel('Count','FontSize',12) xlabel('Month','FontSize',12)%% 年頻數折線圖figureplot(y,D_fre_yr,'k','linewidth',1.5,'Marker','.','MarkerSize',15); set(gca,'FontSize',16) %坐標範圍set(gca,'xminortick','on','yminortick','on');%小刻度打開axis([1949,2018,0,12]); xlabel('Year','Fontsize',18);ylabel('Count','Fontsize',18);grid offbox on%% 距平圖figurebar(y,D_fre_anom,0.8,'FaceColor',[0.4 0.4 0.4]);set(gca,'FontSize',16) %坐標範圍set(gca,'xminortick','on','yminortick','on');%小刻度打開axis([1949,2018,-6,6]); xlabel('Year','Fontsize',18);ylabel('Count','Fontsize',18);grid offbox on%% 挑選強颱風ITC = zeros(size(A));for i=1:length(a1)-1 for j=1:a1(i+1)-a1(i)-1 if A(a1(i)+j,6)==max(A(a1(i)+1:a1(i+1)-1,6))&& A(a1(i)+j,6)>58 ITC(a1(i):a1(i+1)-1,:)=A(a1(i):a1(i+1)-1,:); %強颱風最佳路徑 break end end end%%m=1; for i=1:length(ITC(:,1)) if ITC(i,1)~=0 ITC_bst(m,:) = ITC(i,:); m=m+1; end end %%clear a1 B C lat lon lat1 lon1[a1]=find(ITC_bst(:,1)==66666);B=[];C=[];j=1;l=1;for k=1:length(a1) if 220>ITC_bst(a1(k)+1,3)&&ITC_bst(a1(k)+1,3)>50&&1200>ITC_bst(a1(k)+1,4)&&ITC_bst(a1(k)+1,4)>1050 B(j,:)=ITC_bst(a1(k)+1,:);% 南海生成的強颱風， 未有強颱風在南海生成，所以B為空矩陣 j=j+1; else C(l,:)=ITC_bst(a1(k)+1,:);% 非局地強颱風 l=l+1; endend%% 強颱風生成位置% lat=B(:,3)/10;% lon=B(:,4)/10;lat1=C(:,3)/10;lon1=C(:,4)/10;figurem_proj('Equidistant Cylindrical','lat',[0 40],'long',[100 180]); %投影方式及繪圖範圍設定m_gshhs_l('patch',[0.98 0.98 0.98],'EdgeColor',[0 0 0]); %線條及色塊顏色設定% m_grid(); %邊框及網格設定%m_grid('linestyle','none','fontsize',18);m_grid('box','fancy','tickdir','in');hold on% m_plot(lon,lat,'.k'); m_plot(lon1,lat1,'.r');m_plot(repmat(120,1,18),5:22,'color','m','linewidth',1.5)m_plot(repmat(105,1,18),5:22,'color','m','linewidth',1.5)m_plot(105:5:120,repmat(22,1,4),'color','m','linewidth',1.5)m_plot(105:5:120,repmat(5,1,4),'color','m','linewidth',1.5)2.3 等值線
clc;clear;close alla=peaks(30);[c,h]=contour(a);clabel(c,h,'fontsize',12)
2.4 surf
clc;clear;close allfigure(1)[X,Y] = meshgrid(1:0.5:10,1:20);Z = sin(X) + cos(Y);surf(X,Y,Z)figure(2)[X,Y] = meshgrid(1:0.5:10,1:20);Z = sin(X) + cos(Y);C = X.*Y;surf(X,Y,Z,C)colorbarfigure(3)[X,Y,Z] = peaks(25);CO(:,:,1) = zeros(25); % redCO(:,:,2) = ones(25).*linspace(0.5,0.6,25); % greenCO(:,:,3) = ones(25).*linspace(0,1,25); % bluesurf(X,Y,Z,CO)figure(4)[X,Y] = meshgrid(-5:.5:5);Z = Y.*sin(X) - X.*cos(Y);s = surf(X,Y,Z,'FaceAlpha',0.5);s.EdgeColor = 'none';2.5 contourfclc;clear;close allfigure(1)Z = peaks;contourf(Z)figure(2)x = linspace(-2*pi,2*pi);y = linspace(0,4*pi);[X,Y] = meshgrid(x,y);Z = sin(X) + cos(Y);contourf(X,Y,Z,10)figure(3)[X,Y,Z] = peaks(50);contourf(X,Y,Z,[2 3],'ShowText','on')figure(4)[X,Y,Z] = peaks;contourf(X,Y,Z,[2 2])figure(5)[X,Y,Z] = peaks;contourf(X,Y,Z,'--')figure(6)Z = peaks;[M,c] = contourf(Z);c.LineWidth = 3;figure(7)Z = peaks;Z(:,26) = NaN;contourf(Z)
2.6 pcolorclc;clear;close allfigure(1)X = [1 2 3; 1 2 3; 1 2 3];Y = X';mymap = [1 0 0; 0 1 0; 0 0 1; 1 1 0; 0 0 0];C = [3 4 5; 1 2 5; 5 5 5];pcolor(X,Y,C)colormap(mymap)figure(2)C = hadamard(20);pcolor(C)colormap(gray(2))axis ijaxis squarefigure(3)C = [1 2 3; 4 5 6; 7 8 9];s = pcolor(C);s.EdgeColor = [1 0.7 0.3];s.LineWidth = 6;figure(4)C = [5 13 9 7 12; 11 2 14 8 10; 6 1 3 4 15];s = pcolor(C);s.FaceColor = 'interp';figure(5)[X,Y] = meshgrid(1:20);LY = log(Y);colorscale = [1:20; 20:-1:1];C = repmat(colorscale,10,1);s = pcolor(X,LY,C);tickvals = LY(2:2:20,1)';set(gca,'YTick',tickvals);figure(6)[X,Y] = meshgrid(-3:6/17:3);XX = 2*X.*Y;YY = X.^2 - Y.^2;colorscale = [1:18; 18:-1:1];C = repmat(colorscale,9,1);pcolor(XX,YY,C);figure(7)% 從 R2019b 開始，您可以使用 tiledlayout 和 nexttile 函數顯示分塊繪圖% tiledlayout(1,2)% ax1 = nexttile;C1 = rand(20,10);pcolor(ax1,C1)% ax2 = nexttile;C2 = rand(50,10);pcolor(ax2,C2)
3 Python
3.1 matplotlib+cartopy畫WRF小時降水
from netCDF4 import Datasetimport xarray as xrimport numpy as npimport matplotlib.pyplot as pltimport cartopy.crs as ccrsfrom matplotlib.cm import get_cmapimport matplotlib.colors as colorsfrom cartopy.io.shapereader import Readerimport cartopy.feature as cfeimport cmapsfrom matplotlib.backends.backend_pdf import PdfPagesfrom wrf import (get_cartopy, latlon_coords, to_np, cartopy_xlim, cartopy_ylim,getvar, ALL_TIMES) f = Dataset('wrfout_d01')rainc = getvar(f, 'RAINC', timeidx=ALL_TIMES)rainnc = getvar(f, 'RAINNC', timeidx=ALL_TIMES)tot = rainc.values + rainnc.valueswrf_time = f.variables['Times']dim_tot = tot.shapeh1prep = tot[2:dim_tot[0]:2,:,:] - tot[0:dim_tot[0]-2:2,:,:]lats = f.variables['XLAT'][0]lons = f.variables['XLONG'][0]cart_proj = get_cartopy(rainc)levels = [0.1,0.5,1.,2.,3.,4.,5.,6.,8.,10.,20.,40]norm = colors.BoundaryNorm(boundaries=np.array(levels), ncolors=len(levels)-1)shp_path = '/Users/james/Documents/GitHub/py_china_shp/'reader = Reader(shp_path + 'Province_9/Province_9.shp')provinces = cfe.ShapelyFeature(reader.geometries(), ccrs.PlateCarree(), edgecolor='k', facecolor='none')with PdfPages('multipage_pdf.pdf') as pdf: for i in range(3): fig = plt.figure(figsize=(14,8.5)) ax = fig.add_subplot(111, projection=cart_proj) print(ax) ax.add_feature(provinces, linewidth=0.6) ax.coastlines('50m', linewidth=0.8) pcm = ax.contourf(lons, lats, h1prep[i], levels, norm=norm, transform=ccrs.PlateCarree(), extend='both', cmap=cmaps.cmorph[1:]) pcm.cmap.set_under(cmaps.cmorph.colors[0]) ax.set_title("WRF Prepcipitation {}".format(str(wrf_time[i],'utf-8'))) ax.gridlines(linewidth=0.5, color='gray', linestyle='--') plt.colorbar(pcm, ax=ax, ticks=levels, extendfrac='auto', aspect=18, format='%3.1f', shrink=.95) pdf.savefig() plt.close() 3.2 subplots
import matplotlib.pyplot as pltimport numpy as npt = np.arange(0.0, 2.0, 0.01)s1 = np.sin(2 * np.pi * t)s2 = np.exp(-t)s3 = s1 * s2fig, axs = plt.subplots(3, 1, sharex=True)fig.subplots_adjust(hspace=0)axs[0].plot(t, s1)axs[0].set_yticks(np.arange(-0.9, 1.0, 0.4))axs[0].set_ylim(-1, 1)axs[1].plot(t, s2)axs[1].set_yticks(np.arange(0.1, 1.0, 0.2))axs[1].set_ylim(0, 1)axs[2].plot(t, s3)axs[2].set_yticks(np.arange(-0.9, 1.0, 0.4))axs[2].set_ylim(-1, 1)plt.savefig('F:/Rpython/lp9/plot1.png',dpi=1200)plt.show()import matplotlib.pyplot as pltimport numpy as npfig, axs = plt.subplots(2, 2)cm = ['RdBu_r', 'viridis']for col in range(2): for row in range(2): ax = axs[row, col] pcm = ax.pcolormesh(np.random.random((20, 20)) * (col + 1),cmap=cm[col]) fig.colorbar(pcm, ax=ax)plt.savefig('F:/Rpython/lp9/plot2.png',dpi=1200)plt.show()import matplotlib.pyplot as pltimport numpy as npfig, axs = plt.subplots(2, 2)cm = ['RdBu_r', 'viridis']for col in range(2): for row in range(2): ax = axs[row, col] pcm = ax.pcolormesh(np.random.random((20, 20)) * (col + 1), cmap=cm[col]) fig.colorbar(pcm, ax=axs[:, col], shrink=0.6)plt.savefig('F:/Rpython/lp9/plot3.png',dpi=1200)plt.show()import matplotlib.pyplot as pltimport numpy as npfig, axs = plt.subplots(3, 3, constrained_layout=True)for ax in axs.flat: pcm = ax.pcolormesh(np.random.random((20, 20)))fig.colorbar(pcm, ax=axs[0, :2], shrink=0.6, location='bottom')fig.colorbar(pcm, ax=[axs[0, 2]], location='bottom')fig.colorbar(pcm, ax=axs[1:, :], location='right', shrink=0.6)fig.colorbar(pcm, ax=[axs[2, 1]], location='left')plt.savefig('F:/Rpython/lp9/plot4.png',dpi=1200)plt.show()3.3 Python截取不規則區域內格點
'''author：Masterpiece'''import matplotlib.pyplot as pltfrom mpl_toolkits.basemap import Basemapimport numpy as npimport netCDF4 as ncfrom matplotlib.path import Pathimport shapefileaf=nc.Dataset('F:/Rpython/lp9/grid_within/test.nc')lat=af.variables['latitude'][:]lon=af.variables['longitude'][:]p72=af.variables['p72.162'][0,:,:]/100lat_num=len(lat)lon_num=len(lon)sf=shapefile.Reader("F:/Rpython/lp9/grid_within/Sichuan_province")Shapes = sf.shapes()boundary = Shapes[0].pointsp= Path(boundary)lons, lats = np.meshgrid(lon, lat)grib=np.stack((lons,lats),axis=-1) mask=p.contains_points(grib.reshape(lat_num*lon_num,2)) mask=mask.reshape(lat_num,lon_num)print(lons[mask])lon_leftup=90-2.5/2;lat_leftup=40+2.5/2lon_rightdown=120+2.5/2;lat_rightdown=20-2.5/2fig, ax = plt.subplots(figsize=(14,9))m = Basemap(projection='cyl', llcrnrlat=lat_rightdown, urcrnrlat=lat_leftup, llcrnrlon=lon_leftup, urcrnrlon=lon_rightdown, resolution='i')m.drawcoastlines(linewidth=0.3, color='black')m.readshapefile('F:/Rpython/lp9/grid_within/SC', 'sc', color='k', linewidth=0.6)x,y=m(lons[mask],lats[mask])m.scatter(x,y,s=50, marker='o', color='g', zorder=6)x,y=m(lons[~mask],lats[~mask])m.scatter(x,y,s=50, marker='x', color='r', zorder=6)parallels = np.arange(20,45,5) m.drawparallels(parallels,labels=[True,False,False,False],linewidth=0.6,dashes=[1,4])meridians = np.arange(90,125,5) m.drawmeridians(meridians,labels=[False,False,False,True],linewidth=0.6,dashes=[1,4])plt.yticks(parallels,len(parallels)*[''])plt.xticks(meridians,len(meridians)*[''])plt.savefig('F:/Rpython/lp9/plot5.png',dpi=1200)plt.show()3.4 Python cartopy
'''author：南京大學gavin博士'''import numpy as npimport xarray as xrimport cartopy.crs as ccrsimport cartopy.feature as cfeatfrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERimport matplotlib.pyplot as pltds = xr.open_dataset('F:/Rpython/lp9/EC-Interim_monthly_2018.nc')temp = (ds['t2m'] - 273.15).mean(dim='time') temp.attrs['units'] = 'deg C' proj = ccrs.PlateCarree() fig = plt.figure(figsize=(9,6)) ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1)ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=1) ax.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=1) ax.add_feature(cfeat.LAKES.with_scale('50m'), zorder=1) gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1.2, color='k', alpha=0.5, linestyle='--')gl.xlabels_top = False gl.ylabels_right = False gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER cbar_kwargs = { 'orientation': 'horizontal', 'label': '2m temperature (℃)', 'shrink': 0.8, 'ticks': np.arange(-30,30+5,5)}levels = np.arange(-30,30+1,1)temp.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r', cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())plt.savefig('F:/Rpython/lp9/plot6.png',dpi=1200)plt.show()
3.5 cartopy回歸分析+相關分析'''author：好久不見'''import numpy as np import xarray as xr import os, cmapsfrom sklearn.feature_selection import f_regressionimport matplotlib.pyplot as plt import cartopy.crs as ccrs from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterimport cartopy.io.shapereader as shpreaderfrom cartopy.util import add_cyclic_pointwith xr.open_dataset('/mnt/e/research/data/seasonal/DJF/2.5x2.5/sst.DJF.mean.anom.nc') as f1: pre = f1['sst_anom'][:-1,:,:] lat, lon = f1['lat'], f1['lon']pre2d = np.array(pre).reshape(pre.shape[0],pre.shape[1]*pre.shape[2])del prewith xr.open_dataset('/mnt/e/research/data/seasonal/DJF/2.5x2.5/pc.DJF.sst.nc') as f2: pc = f2['pc'][0,:]del f1, f2A = np.vstack([pc, np.ones(len(pc))]).Tpre_reg = np.linalg.lstsq(A, pre2d)[0][0].reshape(len(lat),len(lon))pre_cor = np.corrcoef(pre2d.T,pc)[:-1,-1].reshape(len(lat),len(lon))pre_cor_sig = f_regression(np.nan_to_num(pre2d), pc)[1].reshape(len(lat),len(lon))area = np.where(pre_cor_sig < 0.01)del pre2d, pc, A, pre_cor_sigpre_reg_cyc, lon_cyc = add_cyclic_point(pre_reg, coord = lon)pre_cor_cyc = add_cyclic_point(pre_cor)nx, ny = np.meshgrid(lon_cyc, lat)del pre_reg, pre_cor, lat, lon, lon_cycplt.figure(figsize = (12, 10))plt.subplots_adjust(hspace = 0.3)ax1 = plt.subplot(211, projection = ccrs.PlateCarree(central_longitude = 180))ax1.coastlines(lw = 0.6)ax1.set_global()c1 = ax1.contourf(nx, ny, pre_reg_cyc, np.arange(-1.5,1.6,0.1), cmap = cmaps.BlWhRe, transform = ccrs.PlateCarree())plt.colorbar(c1, shrink = 1.0, pad = 0.01)sig1 = ax1.scatter(nx[area], ny[area], marker = '.', s = 1, c = 'k', alpha = 0.6, transform = ccrs.PlateCarree())plt.title('Reg.', fontsize = 20)ax1.set_xticks(np.arange(0,361,30), crs = ccrs.PlateCarree())ax1.set_yticks(np.arange(-90,90,15), crs = ccrs.PlateCarree())ax1.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label = False))ax1.yaxis.set_major_formatter(LatitudeFormatter())ax1.set_extent([0,361,-40,85], crs = ccrs.PlateCarree())ax2 = plt.subplot(212, projection = ccrs.PlateCarree(central_longitude = 180))ax2.coastlines(lw = 0.6)ax2.set_global()c2 = ax2.contourf(nx, ny, pre_cor_cyc, np.arange(-1.0,1.1,0.1), cmap = cmaps.BlWhRe, transform = ccrs.PlateCarree())plt.colorbar(c2, shrink = 1.0, pad = 0.01)sig2 = ax2.scatter(nx[area], ny[area], marker = '.', s = 1, c = 'k', alpha = 0.6, transform = ccrs.PlateCarree())plt.title('Cor.', fontsize = 20)ax2.set_xticks(np.arange(0,361,30), crs = ccrs.PlateCarree())ax2.set_yticks(np.arange(-90,90,15), crs = ccrs.PlateCarree())ax2.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label = False))ax2.yaxis.set_major_formatter(LatitudeFormatter())ax2.set_extent([0,361,-40,85], crs = ccrs.PlateCarree())plt.show()
4 徒手轉發DEM連結：https://search.asf.alaska.edu/
SRTM 30m 連結：https://pan.baidu.com/s/1nV0r9RfXvI21v9PraycUYQ 提取碼：7dmv
GDEMV2 連結：https://pan.baidu.com/s/1_G39Rx4j3YVGJH947-NfsA 提取碼：8hp1
GDEMV3 連結：https://pan.baidu.com/s/1UjVuUY6Di3y3SLUUSWvoIg 提取碼：08f5