get.tracks <- function(locations, N=10, x=best.use, se=c(1851, 2010), field="WmaxS", umin=34) { #locations is a data frame with three columns lon, lat, R #R is the search radius in nautical miles (nmi) to look for tracks. #N is the number of tracks to find. #x is the hurricane data file (e.g. best.use). #Each row in the locations database is a point associated with a historical track. #The algorithm searches for all paths that fall within the given radii and generates a score based #on the square root of the sum of the euclidean distances for all storms, from the hourly interpolations. #The output contains of a single list. The first item is a vector of distances, and the second item is a list in which each corresponding #element is a track, for example out$dist[4] is the distance of the 4th possible track, and out$track[[4]] is the actual track from the track file SidDist<-NULL for(i in 1:nrow(locations)) { SidDist<-getSid(lon=locations[i,"lon"],lat=locations[i,"lat"],radius=locations[i,"R"],SidDist=SidDist,se=se,field=field,umin=umin,x=x) } Sid<-SidDist$Sid dist<-apply(SidDist[,-1,drop=F],1,function(x) sqrt(sum(x^2))) od<-order(dist) SidDistO<-data.frame(Sid=Sid,Dist=dist) #SidDistO, not SidDist0 SidDistO<-SidDistO[od,] NN<-nrow(SidDistO) if(NN>=N) SidDistO<-SidDistO[1:N,] else { warning(paste("We found only",NN,"tracks")) N<-NN } tracks<-lapply(SidDistO$Sid,function(i,xx=x) x[x$Sid==i,]) return(list(tracks=tracks,SidDist=SidDistO,N=N,locations=locations)) } getSid<-function(lon,lat,radius,x,SidDist=NULL,se=NULL,field="WmaxS",umin=34) { #x is full dataframe including Sid, if(is.null(SidDist)) { n<-1 } else { x<-x[x$Sid %in% SidDist$Sid,] n<-ncol(SidDist) } xin<-inside.lonlat(x=x,lon=lon,lat=lat,yrlim=se,r=radius) & x[,field] > umin if(sum(xin)==0) stop(paste("Please increase your search radius, no storms selected for",n,"iteration.")) x<-x[xin,] pos<-lonlattoxyz(data=x) pos0<-lonlattoxyz(lon=lon,lat=lat) xdist<-distance.from.center(pos0=pos0,pos=pos) xin<- xdist <= radius xdist<-xdist[xin] if(sum(xin)==0) stop(paste("Please increase your search radius, no storms selected for",n,"iteration.")) x<-x[xin,] dist<-sapply(split(xdist,x$Sid),min) #Get minimum distance Sid<-sapply(split(x$Sid,x$Sid),function(x) x[1]) #Get proper Sid if(is.null(SidDist)) { SidDist<-data.frame(Sid=Sid,dist.1=dist) } else { SidDist<-cbind(SidDist[match(Sid,SidDist$Sid),],dist=dist) colnames(SidDist)[n+1]<-paste("dist",n,sep=".") } SidDist } inside.lonlat <- function(x, lat, lon, yrlim=NULL, r, r0=3440.07) { ###x is a dataframe or matrix with columns "lon" and "lat", ###r is the radius (nmi), r0 is the radius of the earth (nmi) ###returns a subset of xx ###In the merdinal direction r=r0*diff(latlim)/2*pi/180 => difflatlim<-c(-1,1)*r/r0*180/pi latlim<-difflatlim+lat if(is.null(yrlim)) yr<-NULL else yr=x[,"SYear"] ###In zonal direction r=r0*diff(lonlim)/2*pi*cos(lat)/180 ###but cos(lat) changes! since we want r smaller than actual r ###use min(cos(lat)) or min(cos(latlim)) in formula lonlim<-lon+c(-1,1)*r/r0 * 180/pi /min(cos(latlim/180*pi)) xx<-inside.box(lon=x[,"lon"],lat=x[,"lat"],yr=yr,lonlim=lonlim,latlim=latlim,yrlim=yrlim) return(xx) } inside.box<-function(lon,lat,yr=NULL,lonlim,latlim,yrlim=NULL) { if(is.null(yr)){ yr=0 ;yrlim=c(0,0) } ###this function generates a box xlim, ylim ###this fucntion returns a logical vector for each position within the lon,lat, box. return(lon>=lonlim[1] & lon=latlim[1] & lat= yrlim[1] & yr <= yrlim[2] ) } "lonlattoxyz" <- function(lon=NULL, lat=NULL, r = 3440.07,data=NULL ) #3440.07 knots { #This function converts longitude and lattitude 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 earch as 3440.07 nautical miles, #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)) } distance.from.center<-function(pos0,pos) { CX<-pos0[1,1] CY<-pos0[1,2] CZ<-pos0[1,3] X<-pos[,1]-CX Y<-pos[,2]-CY Z<-pos[,3]-CZ return(sqrt(X^2+Y^2+Z^2)) }