栅格地理对象-线性趋势变化图
地理学变量,如降水、气温、污染物浓度等,我们希望计算得到多年栅格图像在空间上的变化趋势、显著水平p、R2、RMSE等等。这里基于matlab和R分别实现线性模型计算的方法和代码。
两者的计算结果一致,效率相当,可任选一种。
Matlab程序
clear; clc; warning off dirpath='K:\liulei\DryNH3\geotif\2surfNH3\yr\'; images=[]; p_intercept=[]; p_slope=[]; fc=[]; fd=[]; fe=[]; p_R2=[]; fg=[]; p_RMSE=[]; p_sig=[]; for i=2008:1:2016 filepath=strcat(dirpath,'\','surfNH3_',num2str(i),'.tif'); img=imread(filepath); [rows, cols]=size(img); temp=reshape(img,rows*cols,1); images=[images temp]; end for i=1:rows*cols i xData=[2008:2016]'; yData= images(i,:)'; if(yData(1)<0) p_intercept=[p_intercept; nan]; p_slope=[p_slope; nan]; p_R2=[p_R2; nan]; p_RMSE=[p_RMSE; nan]; p_sig=[p_sig; nan]; else [m,n]=size(xData); %计算自变量矩阵行列数 X=[ones(m,1),xData]; Y=yData; [B,Bint,E,Eint,Stats]=regress(Y,X); [C P]=corrcoef(xData,yData); %计算相关系数矩阵 R2=C(1,2)^2; [n aa]=size(xData); f=B(1)+B(2)*xData; RMSE=sqrt(sum((f-yData).^2)/n); p_intercept=[p_intercept; B(1)]; p_slope=[p_slope; B(2)]; p_R2=[p_R2; R2]; p_RMSE=[p_RMSE; RMSE]; p_sig=[p_sig; P(1,2)]; end end outpath='K:\liulei\DryNH3\geotif\2surfNH3\yr\Trend\program\matlab\res'; xlswrite([outpath '\parameter.xlsx'],{'p_intercept','p_slope','p_R2','p_RMSE','p_sig'},1,'a1'); xlswrite([outpath '\parameter.xlsx'],p_intercept,1,'a2'); xlswrite([outpath '\parameter.xlsx'],p_slope,1,'b2'); xlswrite([outpath '\parameter.xlsx'],p_R2,1,'c2'); xlswrite([outpath '\parameter.xlsx'],p_RMSE,1,'d2'); xlswrite([outpath '\parameter.xlsx'],p_sig,1,'e2'); a=reshape(p_intercept,rows, cols); b=reshape(p_slope,rows, cols); c=reshape(p_R2,rows, cols); d=reshape(p_RMSE,rows, cols); e=reshape(p_sig,rows, cols); R = georasterref('RasterSize', [rows, cols], ... 'RasterInterpretation', 'cells', ... 'Latlim', [-90 90], 'Lonlim', [-180 180], ... 'ColumnsStartFrom', 'north'); geotiffwrite('K:\liulei\DryNH3\geotif\2surfNH3\yr\Trend\program\matlab\res\p_intercept.tif',a,R); geotiffwrite('K:\liulei\DryNH3\geotif\2surfNH3\yr\Trend\program\matlab\res\p_slope.tif',b,R); geotiffwrite('K:\liulei\DryNH3\geotif\2surfNH3\yr\Trend\program\matlab\res\p_R2.tif',c,R); geotiffwrite('K:\liulei\DryNH3\geotif\2surfNH3\yr\Trend\program\matlab\res\p_RMSE.tif',d,R); geotiffwrite('K:\liulei\DryNH3\geotif\2surfNH3\yr\Trend\program\matlab\res\p_sig.tif',e,R);
R程序
library(maptools) library(ncdf4) library(raster) library(rasterVis) library(RColorBrewer) library(zoo) library('stringr') library('latex2exp') library(marmap) library(sf) library(greenbrown) library('reticulate') ncfname2 <- paste("/home/lei/R/data/GC/surfNH3sate/", paste('surfNH3_', c(2008:2016),".tif", sep="") , sep="") ConvNH4=stack(ncfname2) #ConvNH4=projectRaster(ConvNH4,crs="+proj=longlat +datum=WGS84") ConvNH4v2=as.array(ConvNH4) xyz=dim(ConvNH4v2) ConvNH4v2=array_reshape(ConvNH4v2, c(xyz[1]*xyz[2],xyz[3])) parasFit=c(0) for (index in 1:(xyz[1]*xyz[2])) { index NH3layer=c(0) NH3layer=ConvNH4v2[index,] Height=c(1:9) if(any(is.na(NH3layer))){ parasFit=rbind(parasFit,c(NA,NA,NA,NA)) }else{ fit <- lm(NH3layer~Height) yfitted=predict(fit, data.frame(x=Height)) yfitted=as.numeric(as.character(yfitted)) NH3layerv2=NH3layer R2 <- summary(fit)$r.squared pVal <- anova(fit)$'Pr(>F)'[1] parasFit=rbind(parasFit,c(array(coef(fit)),R2,pVal)) } } #write.csv(parasFit[2:dim(parasFit)[1],],paste('/home/lei/R/data/GC/surface/yr/Guassian',j,'.csv',sep='')) R2layer=array_reshape(parasFit[2:dim(parasFit)[1],3], c(xyz[1],xyz[2])) R2layer <-raster(R2layer) extent(R2layer) <-extent(ConvNH4) crs(R2layer ) <- CRS("+proj=longlat +datum=WGS84") writeRaster(R2layer,paste('/home/lei/R/data/GC/surfNH3sate/Trend/Wet_R2_v3','.tif',sep=""),options=c('TFW=YES'), overwrite=TRUE) print(c(cellStats(R2layer,stat='min'), cellStats(R2layer,stat='max'))) Alayer=array_reshape(parasFit[2:dim(parasFit)[1],1], c(xyz[1],xyz[2])) Alayer <-raster(Alayer) extent(Alayer) <-extent(ConvNH4) crs(Alayer) <- CRS("+proj=longlat +datum=WGS84") writeRaster(Alayer,paste('/home/lei/R/data/GC/surfNH3sate/Trend/Wet_A_v3','.tif',sep=""),options=c('TFW=YES'), overwrite=TRUE) print(c(cellStats(Alayer,stat='min'), cellStats(Alayer,stat='max'))) Blayer=array_reshape(parasFit[2:dim(parasFit)[1],2], c(xyz[1],xyz[2])) Blayer<-raster(Blayer) extent(Blayer) <-extent(ConvNH4) crs(Blayer) <- CRS("+proj=longlat +datum=WGS84") writeRaster(Blayer,paste('/home/lei/R/data/GC/surfNH3sate/Trend/Wet_B_v3','.tif',sep=""),options=c('TFW=YES'), overwrite=TRUE) print(c(cellStats(Blayer,stat='min'), cellStats(Blayer,stat='max'))) player=array_reshape(parasFit[2:dim(parasFit)[1],4], c(xyz[1],xyz[2])) player <-raster(player) extent(player) <-extent(ConvNH4) crs(player) <- CRS("+proj=longlat +datum=WGS84") writeRaster(player,paste('/home/lei/R/data/GC/surfNH3sate/Trend/Wet_p_v3','.tif',sep=""),options=c('TFW=YES'), overwrite=TRUE) print(c(cellStats(player,stat='min'), cellStats(player,stat='max')))