{"word_count": 619, "character_count": 6549}
The National Hurricane Center (NHC) collects datasets with all  storms in North Atlantic, the North Atlantic Hurricane Database (HURDAT). For all sorms, we have the location of the storm, every six jours (at midnight, six a.m., noon and six p.m.). Note that we have also the date, the maximal wind speed - on a 6 hour window - and the pressure in the eye of the storm.It is possible to run the following functionlibrary(XML)extract.track=function(year=2012,p=TRUE){First, we need to extract the list of all stormsloc &#60;- paste("http://weather.unisys.com/hurricane/atlantic/",year,"/index.php",sep="")tabs &#60;- readHTMLTable(htmlParse(loc)) storms &#60;- unlist(strsplit(as.character(tabs[[1]]$Name),split=" "))index &#60;- storms%in%c("Tropical","Storm",paste("Hurricane-",1:6,sep=""),"Depression","Subtropical","Extratropical","Low",paste("Storm-",1:6,sep=""),"Xxx")nstorms  &#60;- storms[!index]and then, for all storms, we go to the appropriate page, and we extract the dataset,TRACK=NULLfor(i in length(nstorms):1){loc=paste("http://weather.unisys.com/hurricane/atlantic/",year,"/",nstorms[i],"/track.dat",sep="")track=read.fwf(loc,skip=3,widths = c(4,6,8,12,4,6,20))names(track)=c("ADV","LAT","LON","TIME","WIND","PR","STAT")track$LAT=as.numeric(as.character(track$LAT))track$LON=as.numeric(as.character(track$LON))track$WIND=as.numeric(as.character(track$WIND))track$PR=as.numeric(as.character(track$PR))track$year=yeartrack$name=nstorms[i]TRACK=rbind(TRACK,track)if(p==TRUE){  cat(year,i,nstorms[i],nrow(track),"\n")}}return(TRACK)}For instance, for 2012, we have&#62; m=extract.track(2012)2012 19 TONY 20 2012 18 SANDY 45 2012 17 RAFAEL 56 2012 16 PATTY 11 2012 15 OSCAR 13 2012 14 NADINE 96 2012 13 MICHAEL 42 2012 12 LESLIE 60 2012 11 KIRK 22 2012 10 JOYCE 12 2012 9 ISAAC 51 2012 8 GORDON 26 2012 7 HELENE 38 2012 6 FLORENCE 21 2012 5 ERNESTO 39 2012 4 DEBBY 18 2012 3 CHRIS 31 2012 2 BERYL 33 2012 1 ALBERTO 20 &#62; head(m)  ADV  LAT   LON      TIME WIND  PR  STAT year name1   1 20.1 -50.8  10/21/18Z   25 1011 LOW 2012 TONY2   2 20.4 -51.2  10/22/00Z   25 1011 LOW 2012 TONY3   3 20.8 -51.5  10/22/06Z   25 1010 LOW 2012 TONY4   4 21.3 -51.7  10/22/12Z   30 1009 LOW 2012 TONYWith a simple loop we can extract all the years (note some manual corrections are necessary, since there are about ten typos in the lists, on the hurricane names).TOTTRACK=NULLfor(y in 2012:1851){TOTTRACK=rbind(TOTTRACK,extract.track(y))}We can visualize all those trajectories,library(maps)map("world",xlim=c(-80,-40),ylim=c(10,50),col="light yellow",fill=TRUE)for(n in unique(TOTTRACK$name)){lines(TOTTRACK$LON[TOTTRACK$name==n],TOTTRACK$LAT[TOTTRACK$name==n],lwd=.5,col="red")}If we look at the distribution of those points, we havelibrary(ks)U=TOTTRACK[,c("LON","LAT")]U=U[!is.na(U$LON),]H=diag(c(.2,.2))fat=kde(U,H,xmin=c(min(U[,1]),min(U[,2])),xmax=c(max(U[,1]),max(U[,2])))z=fat$estimateimage(z)map("world",add=TRUE)Almost a year ago, Christophe Denuse-Baillon - in his actuarial thesis - suggested to use a Markov space process.  I do not want to fit a parametric model, but for all locations, it is possible to use a nonparametric approach. For all locations, on a given grid, it is possible to get all past movements (if any). To be more specific, given a location (a truncated latitude and atruncated longitude), let us extract two informations,the proportion of trajectories which ended at this locationfor all trajectories that did not end at this location, we can extract location 6 hours lated, after being at this location.The code will bemtransition=function(i,j){MOVEMENT=NAsumstop=0p=1idx=which( (TOTTRACK$LON&#62;=gridx[i])&#38;(TOTTRACK$LON&#60;gridx[i+1]) &#38;           (TOTTRACK$LAT&#62;=gridy[j])&#38;(TOTTRACK$LAT&#60;gridy[j+1]) )if(length(idx)&#62;0){MOVEMENT=data.frame(LON=rep(NA,length(idx)),LAT=rep(NA,length(idx)),D=rep(NA,length(idx)))for(s in 1:length(idx)){stops=TRUEif((is.na(TOTTRACK$name[idx[s]+1])==FALSE)&#38;(is.na(TOTTRACK$name[idx[s]])==FALSE)){stops=TOTTRACK$name[idx[s]+1]!=TOTTRACK$name[idx[s]]}if(stops==TRUE ){sumstop=sumstop+1}if(stops==FALSE){d=(TOTTRACK$LON[idx[s]+1]-TOTTRACK$LON[idx[s]])^2+(TOTTRACK$LAT[idx[s]+1]-TOTTRACK$LAT[idx[s]])^2locx=locy=NAif((d&#60;90)&#38; TOTTRACK$LON[idx[s]+1]&#60;0){locx=floor(TOTTRACK$LON[idx[s]+1]/pasgrid)locy=floor(TOTTRACK$LAT[idx[s]+1]/pasgrid)}MOVEMENT[s,]=c(locx,locy,d)}}p=sumstop/length(idx)}return(list(listemouv=MOVEMENT,probstop=p))}For convenience, for all values on a grid, let us store all those values in a big listListTransMat=list()ListStopMat =list()for(i in 1:(length(gridx)-1)){for(j in 1:(length(gridy)-1)){nom=abs(gridx[i])*1000+abs(gridy[j])M=mtransition(i,j)ListTransMat[[nom]] = M$listemouvListStopMat[[nom]]  = M$probstop}}Now we can start some code to generate a trajectory. First, we need a starting point, from one we did observe.n=nrow(TOTTRACK)idx=which(TOTTRACK$name[1:(n-1)]!=TOTTRACK$name[2:n])STARTINGVALUE=TOTTRACK[c(1,idx+1),c("LON","LAT")]If we plot them on the map above, we haveabline(v=gridx,col="white",lwd=.5)abline(h=gridy,col="white",lwd=.5)points(STARTINGVALUE$LON,STARTINGVALUE$LAT,pch=19,col="blue")From that location we look at possible trajectories, and then, we generate two random variables. One to see if we move, or if the trajectory ended where we were. And in the case we move, we drawrandomly among past trajectories. The idea is the following. Start at some specific locationand  then, using the previous function, look were you might be, 6 hours later,Draw one possible location, from those observe previously (note that we might reamin in the same square)sim.mouv=function(LOC){lon=LOC[1]lat=LOC[2]nom=as.numeric(abs(lon)*1000+abs(lat))M=ListTransMat[[nom]]M=M[!is.na(M$LON),]stop=ListStopMat[[nom]]u=runif(1)if(u&#60;=stop) LOC2=NAif(u&#62;stop) LOC2=M[sample(1:nrow(M),size=1),c("LON","LAT")]return(LOC2)}And then, from that new locationwe move again (according to some previous trajectory, observed in the past)Actually, keep in mind that it might also be possible to stop here. The code is - more or less - the followingsim.traj=function(){i=sample(1:nrow(STARTINGVALUE),size=1)LOC=round(STARTINGVALUE[i,])MLOC=LOCstop=FALSEwhile(stop==FALSE){LOC2=sim.mouv(LOC)if(is.na(LOC2))  stop=TRUEif(!is.na(LOC2)) LOC=LOC2; MLOC=rbind(MLOC,LOC)}return(MLOC)}Below, we can visualize two possible trajectories, for hurricanes that we did generate,Nice, isn't it? Now we need to add additional information, about the strength of the hurricaine, again using past experience... But let's keep that for another post!