1

There are tens of ways to cluster spatial points in R but I couldn't find any for spatial lines. What I would like to get is the average trajectory directions that each of these cyclones (spatial lines in my case) has over the territory of France similar to the possibility that trajCluster function from openAir library offers.

library(dplyr) 
library(maptools)
library(geosphere)
library(rgeos)

#load the world boundaries polygon
data(wrld_simpl)
contour <- wrld_simpl[wrld_simpl$NAME =='France',] #only the polygon for France

#download the data
x <- read.table("https://forms.naturwissenschaften.ch/imilast/_ERAinterim_1.5_1979_MTEX/ERAinterim_1.5_NH_M02_19790101_20121231_MTEX.txt?_ga=2.18919096.1825595846.1546710263-1112023567.1546710263", sep="", fill = T, nrows = 20000,
                header = F, skip = 2) # download only the first 20000 rows
colnames(x) <- c("Code","CycloneNo","StepNo","DateI10","Year","Month","Day","Time","LongE","LatN","Intensity1","Intensity2","Intensity3")

x <- x %>% filter(Code!=90) ##Just cleaning to remove those 'NA' rows
x <- x[x$LongE > -120 & x$LongE <= 120,] #subset only the traj between -120 and 120 longE
x$DateI10 <-as.POSIXct(paste(x$Year, x$Month, x$Day, x$Time, sep="-"), "%Y-%m-%d-%H", tz = 'GMT')

# create an empty Spatial Line to be filled later
sline <- SpatialLines(LinesList = list(Lines(Line(matrix(0, ncol = 2)), ID = NA)))
proj4string(sline) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
sline$time <- NA
sline$start <-NA
sline$distance <- NA
sline$mean <- NA

#calculate the lifetime for each track; keep only trajectories for France
for (k in unique(x$CycloneNo)){
  message(which(unique(x$CycloneNo) == k), " from ", length(unique(x$CycloneNo)))
  txt2 <- x[x$CycloneNo==k,]
  if (nrow(txt2) > 1){
    d <- SpatialLines(list(Lines(Line(cbind(txt2$LongE, txt2$LatN)), ID= k)))
    proj4string(d) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
    d$time <- difftime(max(txt2$DateI10), min(txt2$DateI10), units = 'hours')
    d$start <-  as.character(min(txt2$DateI10)) # get the formation date
    d$distance <- lengthLine(d)/1000 # calculate distance
    d$mean <- mean(txt2$Intensity1) # average the sea level surface pressure
    if (gIntersects(d, contour, byid=TRUE) == TRUE){
      sline <- rbind(sline, d) #keep only trajectories crossing the France's territory
    }
  }
}

spplot(sline['time'],sp.layout = list(contour, lwd = 3), 
scales  = list(draw=T), main = 'Mean lifetime (h)') #plot the tracks

### what I would like to get is something like
library(openair)

traj <- importTraj(site = "paris", year = 2009)

## calculate clusters
traj <- trajCluster(traj, method = "Angle", type = "season",n.clusters = 4)

Is it possible to get something like this using my data? is clustering method appropriate here or I should some else (like PCA)?

  • It sounds like you just want to create a grouping column based on a timestamp, correct? This is a far cry from cluster analysis and would be query based and not a statistical model. – Jeffrey Evans Aug 1 at 12:58
  • Well, it would be possible to calculate the mean value for each lat/long point and for each season (using the month column) - resulting one line. But I would need to see the frequency of each pattern so that's why I was thinking to a cluster analysis. – Andrei Niță Aug 1 at 13:06
  • If you are trying to cluster them based on spatial proximity and similarity, you could start by clustering using the Hausdorff distance between all pairs of features. gDistance in package rgeos will compute this for you once you've converted that data frame to a spatial object. – Spacedman Aug 1 at 14:14
  • I wonder if you've tried the new package, 'Rsharpie'? ;-) – Henriette Jager Sep 10 at 1:49
2

Let's work through a Haussdorf clustering of lines.

We'll use the sf package for spatial data and distance calculations:

library(sf)

starting with your final x, lets group everything by cyclone number, make line features, and keep the number of points in the group:

cyclones = x %>% group_by(CycloneNo) %>% mutate(n=n()) %>% summarize(n=mean(n),geometry=st_sfc(st_linestring(cbind(LongE, LatN))))

we keep the number of points because we can't make valid line features with one point and there are some single point cyclones. Let's drop them,

cyclones = cyclones[cyclones$n >1,]

Now make that data frame into a spatial data frame with a lat-long coordinate system:

cyclones = st_sf(cyclones)
st_crs(cyclones)=4326

Make a distance matrix using Haussdorf distance:

dmat = st_distance(cyclones, which="Haussdorf")

And feed that into a clustering function:

cl = hclust(as.dist(dmat))

Use cutree to find the five-cluster point and assign each track with a cluster ID. Plot...

cyclones$class5 = as.factor(cutree(cl, 5))
plot(cyclones[,"class5"])

enter image description here

Which looks interesting. Some cyclones seem to have at least part of themselves in another cluster, but that's the way the Haussdorf distance clustering works.

  • Hi @Spacedman. While this is a great tool for clustering lines, do you know what I could use in order to group the trajectories by their main directions? In the end I need to have several lines (like 3 - 4) and their relative frequency (as in the openair package - there is a function there called trajCluster which groups the HYSPLIT trajectories by their directions). – Andrei Niță Aug 3 at 12:23
  • Its not clear to me yet what exactly you want to do. This answer is clustering lines by a distance measure analogous to clustering points or any other set of variables by a distance measure, and I can't think of any other general principle for clustering things. You need a distance measure and a clustering algorithm. It sounds like you want to cluster by some function of the lines, such as direction. Is the location important for your clusters? ie should all westerly cyclones make a cluster even if they are far apart geographically? If you could edit your Q for clarity that will help. – Spacedman Aug 3 at 12:39
  • thank you Spacedman. I edited the code to be more relevant to my problem. – Andrei Niță Aug 3 at 14:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.