#source("c:/mine/starter.r"); rung("C:/mine/blog/TempLS/complete","LSpic.r") makeSH=function(la,lo,N){#Forms spherical harmonics YLM . # la, lo are lat/lon in degrees n=function(L,M){L*L+2*M+1} # maps from L,M to linear array H=array(0,c(length(la),n(N,N))) # will be SH array p2=pi/180 ss=sin(la*p2); cs=cos(la*p2); #cs>0) G=matrix(0.5/sqrt(pi),length(la),N+1); for(M in 1:N){ # Get Latitude associated Legendres by recurrence G[,M+1]=-cs*(2*M-1)*G[,M] } for(M in 0:N){ g=0; h=H[,n(M,M)]=G[,M+1]; for(L in M:N){ # up diagonal if(L>=N)break; H[,n(L+1,M)]=q=((2*L+1)*ss*h-(L+M)*g)/(L-M+1); #print(c(L,M,q[90],h[90],g[90])) g=h;h=q; } } c0=cos(lo*p2); s0=sin(lo*p2); cc=1;ss=0; for(M in 0:N){ # now longitude trigs for(L in M:N){ J=n(L,M); H[,J]=H[,J]*sqrt((L+L+1)*(1+(M>0))*factorial(L-M)/factorial(L+M)); if(M>0){H[,J-1]= cc*H[,J]; H[,J]= ss*H[,J];} } s=ss; ss=ss*c0+cc*s0; cc=cc*c0-s*s0; # recurrence for trig } H } r=x-rep(L,d[3])-rep(G,each=d[1]) nw=sum(ow) n=1:d[1]+d[1]*(nw-1) o=w[n]>0; no=sum(o) oc=(1:d[1]<7281)[o] N=12 today=Sys.Date() iu=cbind((iv[o,2]+180)%%360-180,iv[o,1]) s1=makeSH(iu[,2],iu[,1],N) ga=rep(-44:44*2,each=180) go=rep(-90:89*2,89) s2=makeSH(ga,go,N) ew=t(s1)%*%s1 y5=r[n][o] w5=w[n][o] y5=y5+G[5,116]-sum(y5*w5)/sum(w5) y1=y5%*%s1 er=s2%*%solve(ew,c(y1)) cl=c("#7700ff","#4094ff","#78ccff","#99eefe","#dcffdc", "#ffffff","#ffff00","#ffcc00","#ff7f00","#ff0000","#7f0000"); lv=c(-4,-2,-1,-.5,-.2,.2,.5,1,2,4); lk=outer(er,lv,">") lk=rowSums(lk)+1 require("maps") mp=map(type="n") ik=-5.5:5*32.8 il=rep(-100,11) graphics.off() fs=c("map.png","stats.png") j=-5:0+nw md=month.abb[(j-1)%%12+1] df=paste(md[6],d[3]-101,"d",sep="") png(fs[1],width=800) # Temp map plot(range(go),c(-110,90),type="n",axes=F,xlab="",ylab="",asp=1, main=pc("GHCNV3+ERSST TempLSV3 ",df, " anomalies base 1961-90, made ",today)) rect(go,ga-1,go+2,ga+1,col=cl[lk],border=NA) lines(mp$x,mp$y) rect(ik,il,ik+33,il+4,col=cl) text(ik[-1],il[-1]-10,paste(lv)) dev.off() png(fs[2],width=800) # stations map plot(range(go),range(ga),type="n",axes=F,xlab="",ylab="", main=pc("GHCNV3+ERSSTv4 ",no," stations reporting as at ",today)) points(iu,pch=19,cex=0.5,col=c("red","green")[2-oc]) lines(mp$x,mp$y) dev.off() for(i in 1:2)file.copy(fs[i],paste("O:/data/freq/",df,sep=""),over=T) b=readLines("reptemplate.html") cx=c("@date","@stat","@m0","@t0","@job") cz=c(paste(today),no,pc(md,s="