require(splines) nsd<-function(knotTimes, predTimes, deriv = 1, order = 3, periodic=F) { ###This function returns an interpolation matrix for values at knotTimes ###to values at the predicted times, by solving for the basis coefficients ###and summing the basis function evaluated at the predTimes ###deriv=0, 1 or 2. order not yet implemented lknot <- length(knotTimes) knots <- knotTimes[c(-1, - lknot)] Boundary.knots <- knotTimes[c(1, lknot)] return(ns1(x = predTimes, knots = knots, Boundary.knots = Boundary.knots, intercept = T, deriv = deriv, periodic=periodic) %*% solve(ns1(x = knotTimes, knots = knots, Boundary.knots = Boundary.knots, intercept = T , periodic=periodic))) } ns1<-function(x, df = NULL, knots = NULL, intercept = F, periodic = F, Boundary.knots = range(x), derivs = rep(0, length(x))) { ##This is an updated Splus function to add the "periodic" argument nx <- names(x) x <- as.vector(x) nax <- is.na(x) if(nas <- any(nax)) x <- x[!nax] if(!missing(Boundary.knots)) { if(length(Boundary.knots) == 1 && Boundary.knots == T) { Boundary.knots <- range(knots) knots <- sort(knots)[ - c(1, length(knots))] } Boundary.knots <- range(Boundary.knots) outside <- (ol <- x < Boundary.knots[1]) | (or <- x > Boundary.knots[2]) } else outside <- rep(F, length = length(x)) if(!missing(df) && missing(knots)) { # df = number(interior knots) + 1 + intercept nIknots <- df - 1 - intercept if(nIknots < 0) { nIknots <- 0 warning(paste("df was too small; have used ", 1 + intercept)) } if(nIknots > 0) { knots <- seq(from = 0, to = 1, length = nIknots + 2)[ - c(1, nIknots + 2)] knots <- quantile(x[!outside], knots) } else knots <- NULL } Aknots <- sort(c(rep(Boundary.knots, 4), knots)) if(!missing(derivs) && length(derivs) == 1) derivs <- rep(derivs, length(x)) if(any(outside)) { basis <- array(0, c(length(x), length(knots) + 4)) if(any(ol)) { k.pivot <- Boundary.knots[1] xl <- cbind(1, x[ol] - k.pivot) tt <- spline.des(Aknots, rep(k.pivot, 2), 4, c(0, 1))$design basis[ol, ] <- xl %*% tt if(!missing(derivs) && any(derivs[ol] > 0)) { basis[ol & derivs == 1, ] <- rep(tt[2, ], each = sum(ol & derivs == 1)) basis[ol & derivs > 1, ] <- 0 } } if(any(or)) { k.pivot <- Boundary.knots[2] xr <- cbind(1, x[or] - k.pivot) tt <- spline.des(Aknots, rep(k.pivot, 2), 4, c(0, 1))$design basis[or, ] <- xr %*% tt if(!missing(derivs) && any(derivs[or] > 0)) { basis[or & derivs == 1, ] <- rep(tt[2, ], each = sum(or & derivs == 1)) basis[or & derivs > 1, ] <- 0 } } if(any(inside <- !outside)) basis[inside, ] <- spline.des(Aknots, x[inside], 4, derivs[inside])$design } else basis <- spline.des(Aknots, x, 4, derivs)$design if(periodic) { const12 <- spline.des(Aknots, rep(Boundary.knots, each = 2), 4, c(1, 2, 1, 2))$design const <- const12[3:4, ] - const12[1:2, ] } else const <- spline.des(Aknots, Boundary.knots, 4, c(2, 2))$design if(!intercept) { const <- const[, -1, drop = F] basis <- basis[, -1, drop = F] } qr.const <- qr(t(const)) basis <- as.matrix((t(qr.qty(qr.const, t(basis))))[, - (1:2)]) n.col <- ncol(basis) if(nas) { nmat <- matrix(NA, length(nax), n.col) nmat[!nax, ] <- basis basis <- nmat } dimnames(basis) <- list(nx, 1:n.col) a <- list(degree = 3, knots = knots, Boundary.knots = Boundary.knots, intercept = intercept, class = c("ns", "basis")) attributes(basis) <- c(attributes(basis), a) basis } "xyztolonlatheight"<-function(x=NULL,y=NULL,z=NULL,data=NULL,r0=3440.07) #3440.07 nautical miles, change for other units { #This function is used many places to convert xyz coordintes to lonlatheight, #inputs: x,y,z equal length vectors of x,y,z coordinates, if data argument is given will use data$x, data$y, data$z #r0 is the radius of the earth, here given as 3440.07 nm (mean). #returns matrix with lon, lat, and height columns if(!is.null(data)) { x<-data[,"x",drop=T] y<-data[,"y",drop=T] z<-data[,"z",drop=T] } R<-sqrt(x^2+y^2+z^2) degtorad<-180/pi lat<-degtorad*asin(z/R) lon<-degtorad*atan2(y,x) height<-R-r0 return(cbind(lon=lon,lat=lat,height=height)) } "lonlattoxyz" <- function(lon=NULL, lat=NULL, r = 3440.07, data=NULL) #3440.07 nautical miles { #This function converts longitude and latitude to x,y,z coordiantes. Height is zero. #Input: lon, lat, vectors or data, data frame with lon lat columns, # r is the radius of earth as 3440.07 kt, #returns matrix with columns x,y,z if(!is.null(data)) { lon<-data[,"lon",drop=T] lat<-data[,"lat",drop=T] } rbase <- r * cos((lat * pi)/180) cbind(x = rbase * cos((lon * pi)/180), y = rbase * sin((lon * pi)/180), z = r * sin((lat * pi)/180)) } "xyztouvr"<-function(radii,...) { ### compute u,v velocity vectors ### ... are vector matrices associated with the radii matrix ### radii is an n by 3 matrix of positions associated with list ### returns transformed matrices associated with the radii matrix, ### with column names u,v,w associated with the vertical east and north ### vector components. if(!is.matrix(radii)) { radii<-as.matrix(radii) colnames(radii)<-c("x","y","z") } dv<-dim(radii) if(dv[2]!=3) stop("Invalid positon matrix, must have at least three columns") vecmat<-list(...) w<-names(vecmat) vecmat<-lapply(vecmat,function(x) { if(is.matrix(x)) return(x) else { xx<-as.matrix(x) colnames(xx)<-c("x","y","z") rownames(xx)<-rownames(x) xx } } ) x<-radii[,"x"] y<-radii[,"y"] z<-radii[,"z"] radxy2<-x^2+y^2 radxy<-sqrt(radxy2) R2<-radxy2+radii[,"z"]^2 R<-sqrt(R2) ### The transformations are determined from two rotational forms ### w= 1 [ r 0 z] [ x y 0] Vx ### u= --- [ 0 R 0] [ -y x 0] Vy ### v= R*r [-z 0 r] [ 0 0 r] Vz ### the first matrix rotates the coordinates around the z axis, the second matrix ### rotates the coordinates around the rotated y axis resulting in w the vertical velocity ### u the westerly velocity (eastward moving) and v the southerly velocity (northward moving) wmat<-cbind(x,y,z)/R umat<-cbind(-y,x,0)/radxy vmat<-cbind(-z*x,-z*y,radxy2)/(radxy*R) ruv<-lapply(vecmat,function(v,wmat,umat,vmat) { data=cbind(u=rowSums(umat*v),v=rowSums(v*vmat),w=rowSums(v*wmat)) data=cbind(data,maguv=sqrt(rowSums(data[,c("u","v")]^2)),diruv = 180 * (( atan2(data[,"u"],data[,"v"])/pi) %% 2)) data } ,wmat=wmat,umat=umat,vmat=vmat ) names(ruv) <- names(vecmat) return(ruv) } "interpolate"<-function(split.tracks,deriv=0,deltat=1,order=3,timefield="Shour",progressbar=T) { if(!progressbar || !require(tcltk)) progressbar<=F ###This function generates the information needed for a given derivative of the split tracks ###returns a list of split tracks. ###split.tracks, list of tracks each containing four columns, x,y,z and Sdays. ###return list of x,y,z and Sdays. or the derivatives of x,y,z with one position on each ###hour of the track. Note storms that have missing information still get interpolated on every ###hour, but the velocities should be small. (Multiple event storms, that is) ###order is not used at this time, needs to be incorporated, means a change in the trackinfo R<-length(split.tracks) if(progressbar) pb<-tkProgressBar(title = "Interpolation", label = "Initializing", min = 0, max = R, initial = 0, width = 300) interpolation<-vector(mode="list",length=R) names(interpolation)<-names(split.tracks) interp.fun<-function(x,deltat,deriv,order,timefield) { nrx<-nrow(x) x<-as.matrix(x) rownames(x)<-NULL ## print(nrx) x<-as.matrix(x) timeid<-which(colnames(x)==timefield) if(length(timeid)!=1) stop("Invalid or Missing timefield") times<-x[,timeid] times.interp<-deltat*round(seq(min(times),max(times),deltat)/deltat) #conform to deltat boundaries Need to test... if(nrow(x)==1) { if(deriv==0) return(x) else return( cbind(x[,-timeid,drop=F]*0,x[,timeid,drop=F])) } vmat<-nsd(times,times.interp,deriv=deriv,order=order) interp.data<-vmat%*%x colnames(interp.data)<-colnames(x) ###add better ones later for rowIds(interp.data) interp.data[,timeid]<-times.interp while(F) print(interp.data) return(interp.data) } for(i in 1:R) { interpolation[[i]]<-interp.fun(x=split.tracks[[i]],deltat=deltat,deriv=deriv,order=order,timefield=timefield) if(progressbar) setTkProgressBar(pb,value=i,label=paste(i,"of", R, "storms complete",sep=" ")) } return(interpolation) } "interpolate.tracks"<-function(tracks,hours=c(1),order=3,deriv=1,get.land=F,get.uvw=T,lgrid=landGrid, col.ignore=c("lat","lon","M","hr","Shour","ISOtime"),col.dup=c("Sn","SYear","Yr","Mo","Da","count","Type","name","IBtracSid","Basin","SubBasin"),default.na=list(Pmin=1020,PminS=1020),createindex=F) #This function is run once on the tracks #The function will interpolate all fields unless ignored, or duplicated #a duplicated field is one in which we repeat the previous value between interpolations such that the field value is changed #only if the corresponding interpolated time is greator than or equal to the actual time value of an observation, thus #if we use hours=1, and we have observations every 6 hours we would repeat each row 6 times. #The hr field is special and will be included if it exists as (Shour+hr[1])%%24 no other clock arithmetic will be performed. ###pmin since we are not using these values for track relocation, they #are to be added later, However Wmax will be used since it is widely #available and required for matching tracks. We use all the tracks #in the dataset to come up with a set of interpolated tracks and #track velocities. #This function should be updated for smoothness on tracks. #Inputs: imported Best Track data set (using import.tracks or historical #storm vitals with both days and Sid fields #Outputs: New Best Track Data set with hourly interpreted data and #x,y,z coordinate information to be used with trackinterp or historical storm vitals ####Note there will be no NAs in the output, a column of all zeros for a given track is an indication that ###the value does not exist or if r34=0 or Wmax<34, then r34 does not exist. Similiarly for r50 ###Initially we interpolate all the data, however, for the tracks used, we do not interpret all the data. ###This function returns a list ###consisting of the following components. ###Data for all the tracks with interpolated data ###index, a named list by Storm id. Each list element consists of the rows in Data corresponding to the StormId ###start location, a vector consisting of the start position within Data, so that Data[start,] is start information. ###with NAs for storms in which the start location cannot be identified. ###interpolated lat,lon, x,y,z ###vx,vy,vz; (for deriv=1) (actual velocity is vx,vy, vz + -omega*y,omega*x,0 ###ax, ay, az; (for deriv=2) ###(actual acceleration is ax,ay,az+2*omega*(-vy,vx,0)-omega^2*(x,y,0) rel+coreolis+centrip ### u,v,r components may be calculated from xzytouvr if needed ###hours=1, resolution of model, number of hours between interpolated points ###order, order of interpolation, currently fixed at three, to be updated to five or seven. ###deriv=1, calulate interpolated positions, and velocities only, use xyz coordinates. ###Note, we assume that only the x,y,z and hours components needed for derivatives, { if(hours != 1 && hours != .1) stop("Interpolation only allowd for 1 and .1 hours at this time") fields<-colnames(tracks) tracks.pos<-as.integer(round(tracks$Shour/hours)) col.ignore<-is.element(fields,col.ignore) col.dup<-is.element(fields,col.dup) if(!any(col.dup)) stop("Must duplicate at least one column") col.interp<-!(col.ignore | col.dup) z<-lonlattoxyz(tracks$lon,tracks$lat) rownames(z)<-rownames(tracks) if((ldna<-length(default.na)) > 0) { if(!is.list(default.na) || is.null(names(default.na))) {stop("default.na must be a named list")} default.na.names<-names(default.na) for(i in 1:ldna) { col.default<-is.element(fields,default.na.names[i]) & col.interp if(any(col.default)) tracks[,col.default][is.na(tracks[,col.default])]<-default.na[[i]]+1i } } tracks[,col.interp][is.na(tracks[,col.interp])]<-0+1i imaginary<-sapply(tracks,mode)=="complex" if(any(imaginary)) { imnames<-fields[imaginary] Nadata<-apply(tracks[,imaginary],2,Im) colnames(Nadata)<-paste(imnames,".im",sep="") tracks[,imaginary]<-apply(tracks[,imaginary],2,Re) tracks<-cbind(tracks,Nadata) col.interp<-c(col.interp,rep(TRUE,sum(imaginary))) } ztracks<-split(data.frame(Shour=tracks$Shour,z,tracks[,col.interp,drop=F]),tracks$Sid) ##xdot, ydot and zdot are column names, to differentiate between x, y and z. zvel<-lapply(ztracks,function(z) data.frame(xdot=z$x,ydot=z$y,zdot=z$z,Shour=z$Shour)) zpos<-interpolate(ztracks,deltat=hours,order=order,deriv=0,timefield="Shour") zvel<-interpolate(zvel,deltat=hours,order=order,deriv=1,timefield="Shour") interp.pos<-as.integer(round(unlist(lapply(zpos,"[",,"Shour"))/hours)) interp.Shour<-interp.pos*hours interp.Sid<-as.integer(round(unlist(lapply(zpos,"[",,"Sid")))) Smax<-as.integer(round(10^((trunc(log10(max(tracks.pos))))+1))) #Get easy to use value interp.id<-interp.Sid*Smax+interp.pos tracks.id<-tracks$Sid*Smax+tracks.pos row.map<-simple.map(interp.id,tracks.id) tracks.dup<-tracks[row.map,fields[col.dup],drop=F] tracks.hrstart<-unlist(lapply(split(tracks$hr,tracks$Sid),function(x) x[rep(1,length(x))])) tracks.shourstart<-unlist(lapply(split(tracks$Shour,tracks$Sid),function(x) x[rep(1,length(x))])) interp.hr<-(tracks.hrstart[row.map]+(interp.Shour-tracks.shourstart[row.map]))%%24 interp.zpos<-unsplitS(zpos,addNAs=F) rm(zpos) interp.zvel<-unsplitS(zvel,addNAs=F)[,-4] rm(zvel) rm(ztracks) interp.zpos<-data.frame(interp.zpos) interp.zpos$Sid<-interp.Sid rm(interp.Sid) #print(zpos) lonlatzpos<-xyztolonlatheight(data=interp.zpos) trackinterp<-data.frame(tracks.dup,hr=interp.hr,lonlatzpos, interp.zpos,interp.zvel) trackinterp$ISOtime<-ISOdate(trackinterp$Yr,trackinterp$Mo,trackinterp$Da,trackinterp$hr) ## names(trackinterp)<-as.character(tracks$Sid) ## cbind(best, M = bestM, lonlattoxyz(data = best[, c("lon", "lat")]), WmaxMS = best$Wmax * 0.5144) if(get.land) trackinterp$M<-is.land(grid=lgrid,points=trackinterp[,c("lon","lat")],method="l") if(get.uvw) trackinterp<-cbind(trackinterp,xyztouvr(trackinterp[,c("x","y","z")],trackinterp[,c("xdot","ydot","zdot")])[[1]]) if(createindex) index<-split(1:nrow(trackinterp),trackinterp$Sid) else return(trackinterp) return(list(data=trackinterp,index=index)) } "simple.map"<-function(x,y){ cut(x,c(y,Inf),labels=F,right=F) } unsplitS <-function(x, numeric = T, addNAs = F){ addNAs <- as.numeric(addNAs) ##Add 0,1 or n rows of NAs if(!is.list(x)) stop("Must have list") find.data <- function(x) { length(x) && nrow(x) } use.all <- sapply(x, find.data) x <- x[use.all] dims <- sapply(x, function(x) dim(x)) if(any(dims[2, ] != dims[2, 1])) stop("Invalid Split different number of columns for some splits") rows <- cumsum(dims[1, ] + addNAs) lrows <- length(rows) out <- x[[1]][rep(1, rows[lrows]), , drop = F] out[, ] <- NA rownames(out) <- unlist(lapply(x, function(x, n) c(rep(NA, n), rownames(x)), n = addNAs)) start <- addNAs + 1 + c(0, rows[ - lrows]) end <- rows for(i in 1:lrows) out[start[i]:end[i], ] <- x[[i]] return(out[(start[1]:end[lrows]), ]) } "is.land"<-function(grid=landGrid,points,center=NULL,method="l") { ###This function determines whether a vector of points is inside hull ###by considering the number of rotations a point on the curve makes about a single point. ###The result counts is this number it should be the number of counter clockwise rotations. ###The sign denotes the curve direction (+ couterclockwise) about an interior point. ###This function was converted from a C program ###The center is chosen so that the grid comprises a simple closed curve ###The result is a vector with the number of track crossings ###This works for all closed curves ###Method=1 use lapply, slower but works on large data sets ###Method=2 use outer3, faster, but works only on single tracks. if(!is.null(center)) grid<-rbind(center,grid) xp<-grid[,1] yp<-grid[,2] x<-points[,1] y<-points[,2] if(any(colnames(points)[1:2] != colnames(grid)[1:2])) warning(paste("points column names:",colnames(points)[1:2],"are not equal to","grid column names:",colnames(grid)[1:2])) lx<-length(xp) ##Do not add last point, remove first point instead assume that each loop is closed by having the first and last points ##duplicated. (use close.curves) xpi<-xp[-1] xpj<-xp[-lx] ypi<-yp[-1] ypj<-yp[-lx] gridcross<-cbind(xpi,xpj,ypi,ypj) xy<-cbind(x,y) ##get x crossings if(method=="l") { crossings<-unlist(lapply(1:nrow(xy), function(i,xy,grid) { x<-xy[i,1] y<-xy[i,2] xpi<-grid[,1] xpj<-grid[,2] xi<-x < xpi xj<-x < xpj ### xcross=(xj&(!xi)) | (xi&(!xj) ) == vi<-(xj-xi) vi[is.na(vi)]<-0 ###Closed Curves separated by NAs vil<-as.logical(vi) if(!any(vil)) return(0) vi<-vi[vil] grid<-grid[vil,] xpi<-grid[,1] xpj<-grid[,2] ypi<-grid[,3] ypj<-grid[,4] xdelta<-xpi-xpj ydelta<-ypi-ypj sum(vi*(y<(ypi+(x-xpi)*ydelta/xdelta))) },xy=xy,grid=gridcross)) is.inside<-as.logical(crossings) } else {##Keep for later, may be faster ydelta<-ypj-ypi xdelta<-xpj-xpi gridcross<-cbind(xdelta,ydelta,xpi,ypi) xj<-outer(x,xpj,"<") xi<-outer(x,xpi,"<") ## xcross=(xj&(!xi)) | (xi&(!xj) ) == xcross<-(xj-xi) xcross[is.na(xcross)]<-0 ###Closed Curves separated by NAs crossings<-outer3(x=xy,y=gridcross,fun=function(xy,grid,v) { v<-as.vector(v) x<-xy[,1] y<-xy[,2] xdelta<-grid[,1] ydelta<-grid[,2] xp<-grid[,3] yp<-grid[,4] ifelse(v,v*(y<(yp+(x-xp)*ydelta/xdelta)),0) },v=xcross ) is.inside<-as.logical(rowSums(crossings)) } # print(all(is.insideo==is.inside)) return(is.inside) } ##C code: ## int pnpoly(int npol, float *xp, float *yp, float x, float y) ## { ## int i, j, c = 0; ##for (i = 0, j = npol-1; i < npol; j = i++) { ## if ((((xp[i]<=x) && (x