Acquisizione dei dati

Analizziamo gli ultimi dati relativi alle posizioni delle navi in Adriatico dal sito del progetto SHAPE.

dat <- get_ships()
## OGR data source with driver: GML 
## Source: "/tmp/RtmpDCo8Tu/filebecc3a497f5a", layer: "bs_v_ship_positions_pt"
## with 740 features
## It has 10 fields
## INFO [2018-02-25 15:56:02] Downloaded layer bs_v_ship_positions_pt from http://atlas.shape-ipaproject.eu/wfs
colnames(dat)
##  [1] "x"        "y"        "gml_id"   "gid"      "mmsi"     "velocity"
##  [7] "angle_ms" "date"     "name"     "style"    "type"     "country"

Guardiamo i dati in una tabella interattiva, selezionando solo alcuni campi più interessanti.

library(dplyr)
DT::datatable(dat %>% select(x,y,mmsi,velocity,angle_ms,name,country,type,date),
              filter="top")

Possiamo anche metterli su una mappa. Cliccando sui punti vedremo le informazioni sulle navi.

plot_last_ships(dat)

Controllo dei dati

Possiamo fidarci delle informazioni relative a ciascuna nave? Facciamo qualche controllo incrociato. Per esempio, possiamo calcolare la velocità effettiva e confrontarla con quella dichiarata. Inoltre, i codici MMSI con meno di 9 caratteri e velocità zero sono probabilmente fari o altre infrastrutture costiere. Sarà interessante vedere come erano classificati nel dataset originario. Infine, dalle prime tre cifre dei codici MMSI si può identificare lo Stato di origine di una nave. Carichiamo dunque i dati di alcune giornate (ci vorrà un po’ di tempo) e analizziamoli.

load_ship_data(from = as.Date("2017-12-06"),
               to = as.Date("2018-02-21")) -> dat
ship_speed(dat) -> dat
dat %>%
  group_by(mmsi) %>% 
  mutate(vmean=mean(vel_kmh,na.rm=T)) %>% 
  ungroup() %>%
  mutate(type.new = ifelse(nchar(mmsi)<9 & vmean<0.01,"Coast station",type),
         country.new = mmsi2country(mmsi)) -> dat

Confrontiamo il nuovo campo country.new con quello originale, country.

library(ggplot2)
ggplot(dat %>% group_by(country, country.new) %>% summarize(n=n()),
       aes(x = country, y = country.new, fill = n)) +
  geom_raster() +
  theme_bw() +
  scale_x_discrete(drop = FALSE) +
  scale_y_discrete(drop = FALSE) +
  theme(axis.text.x = element_text(angle = 270, hjust = 0),
        aspect.ratio = 1,
        legend.position = "none")

Si vede che dal codice MMSI abbiamo recuperato più informazioni, senza però contraddire quelle che già avevamo. Ora confrontiamo il nuovo campo type.new con quelli originale, type.

library(ggplot2)
ggplot(dat %>% group_by(type, type.new) %>% summarize(n=n()),
       aes(x = type, y = type.new, fill = n)) +
  geom_raster() +
  theme_bw() +
  scale_x_discrete(drop = FALSE) +
  scale_y_discrete(drop = FALSE) +
  theme(axis.text.x = element_text(angle = 270, hjust = 0),
        aspect.ratio = 1,
        legend.position = "none")

Si vede che le infrastrutture costiere erano in origine classificate “Unknown”. Vediamo ora se la velocità effettiva differisce da quella dichiarata (la convertiamo da nodi a kilometri all’ora).

dat %>% mutate(country=country.new, type=type.new) -> dat
ggplot(dat %>% mutate(vel_orig=as.numeric(velocity)*1.852),
       aes(x=vel_orig, y=vel_kmh)) +
  geom_bin2d(bins=100) +
  scale_x_log10(limits=c(0.1,100), breaks=c(0.2,0.5,1,2,5,10,20,50,100)) + 
  scale_y_log10(limits=c(0.1,100), breaks=c(0.2,0.5,1,2,5,10,20,50,100)) +
  geom_abline(slope=1,intercept=0,linetype=3, col="grey40") +
  scale_fill_gradientn(colours = c("grey70","navy","red")) +
  theme_bw()

Per la maggior parte dei casi c’è una buona corrispondenza, ma non sempre. D’ora in avanti useremo la velocità calcolata.

Conteggio delle navi

Escludiamo i fari e altre infrastrutture costiere, poi contiamo le navi per ciascun paese e per ciascun tipo.

dat %>% filter(type!="Coast station") -> dat
dat %>% group_by(mmsi) %>% select(mmsi,type,country) %>%
  slice(1) %>% ungroup() %>%
  group_by(type,country) %>%
  summarise(n=n()) %>% arrange(desc(n)) %>%
  DT::datatable(filter="top",
                options = list(
                  searching = T,
                  pageLength = 15))

Individuazione dei porti

Di ciascuna nave individuiamo le soste, cioè le fasi in cui la velocità scende sotto 0.5km/h.

library(sp)
dat %>% mutate(still=vel_kmh<0.5) %>% arrange(time) %>%
  group_by(mmsi) %>% 
  mutate(was.still = lag(still),
         change.status = still!=was.still,
         id.status = cumsum(ifelse(is.na(change.status),F,change.status))) %>%
  group_by(mmsi,id.status) %>% 
  summarise(dist_km=round(sum(dist_km)),
            duration = as.numeric(difftime(last(time),first(time), units = "mins")),
            x=median(x),y=median(y),
            vel_kmh=round(dist_km/(duration/60))) %>%
  filter(duration>10) %>% 
  ungroup() %>%
  group_by(mmsi) %>% 
  mutate(status=if_else(vel_kmh==0,"stop","go"),
         change.status = status!=lag(status),
         id.status = cumsum(ifelse(is.na(change.status),F,change.status))) %>%
  group_by(mmsi,id.status) %>%
  summarise(dist_km=sum(dist_km,na.rm=T),
            duration=sum(duration,na.rm=T),
            x=median(x),y=median(y),
            vel_kmh=round(dist_km/(duration/60))) -> dat_aggr
dat_aggr %>% filter(vel_kmh==0) -> dat_pause

Individuiamo le soste di durata superiore a 12 ore. Avverranno probabilmente in corrispondenza di punti di ancoraggio o di approdo.

dat_pause %>% filter(duration>12*60) -> dat_longpause
dat_longpause %>% ungroup() %>% select(x,y) %>% distinct() -> dd

Aggreghiamo questi punti di sosta in 60 cluster, che corrisponderanno a porti o zone di ancoraggio.

library(cluster)
ncl <- 60
kmeans(dd,ncl, nstart=10) -> cc
dat_cluster <- data.frame(dd, sprintf(as.numeric(rownames(fitted(cc))),fmt="%03i"), fitted(cc))
colnames(dat_cluster) <- c("x","y","cluster","x.centroid","y.centroid")

Vediamo la mappa completa dei porti.

library(ggplot2)
MapData <- map_data(map = "world",
                    region = c("Italy","Slovenia","Croatia","Bosnia and Herzegovina"))
sh <- rep(1:5,ncl/5)
names(sh) <- sort(unique(dat_cluster$cluster))
co <- rep(RColorBrewer::brewer.pal(n=5,"Set1"),each=ncl/5)
names(co) <- sort(unique(dat_cluster$cluster))
ggplot() +
  geom_point(data=dat_cluster, aes(x=x.centroid,y=y.centroid,col=cluster,shape=cluster)) +
  scale_x_continuous(limits=c(12,18))+scale_y_continuous(limits=c(40,46)) +
  scale_shape_manual(values = sh)+
  geom_path(data=MapData, aes(x=long, y=lat, group = group),colour="darkgrey" ) +
  theme_bw() -> p
p

Vediamo in dettaglio la zona settentrionale, visualizzando tutti i punti di sosta lunga, non solo i centroidi.

ggplot() +
  geom_point(data=dat_cluster, aes(x=x,y=y,col=cluster,shape=cluster)) +
  scale_x_continuous(limits=c(12,14))+scale_y_continuous(limits=c(45,46)) +
  scale_shape_manual(values = sh)+
  geom_path(data=MapData, aes(x=long, y=lat, group = group),colour="darkgrey" ) +
  theme_bw() + theme(legend.position = "none") -> p
p

La funzione get_geonames, che si appoggia al pacchetto googleway, ci permette di recuperare, per alcuni di questi porti, il nome della località più vicina. Poiché il servizio ha un limite giornaliero di interrogazioni, in fase di sviluppo conviene sospenderlo.

if(exists("ports") && any(ports$name!=ports$cluster)) {
  activ_google <- FALSE
} else {
  activ_google <- TRUE
}
ports <- dat_cluster %>% select(cluster, x.centroid, y.centroid) %>% distinct()
if(activ_google) {
  ports$name <- get_geonames(lon=ports$x, lat=ports$y)
  ports$name[ports$name==""] <- ports$cluster[ports$name==""]
} else {
  ports$name <- sprintf(ports$cluster,fmt = "%03i")
}

Ora la tabella ports contiene nomi e codici dei porti e le coordinate dei centroidi. La colleghiamo alla tabella dei punti di sosta lunga dat_cluster.

dat_cluster <- left_join(dat_cluster %>% 
                           select(x,y,cluster,x.centroid,y.centroid) %>% 
                           distinct() %>% 
                           mutate(cluster=as.character(cluster)),
                         ports %>% 
                           select(cluster,name)%>% 
                           mutate(cluster=as.character(cluster)),
                     all.x=T)

Dai punti di sosta costruiamo i poligoni che li contengono.

ports_polygons <- points2SpatialPolygonsDataFrame(dat_cluster)

Nella mappa interattiva passando il mouse sopra i poligoni possiamo vedere il codice del porto.

library(leaflet)

pal <- colorNumeric(rainbow(ncl), domain = c(1:ncl))
leaflet(data = ports_polygons) %>% addTiles() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  #addCircleMarkers(~x, ~y, fillOpacity = 1,
  addPolygons(color = ~pal(as.numeric(cluster)), 
              stroke = TRUE,
              label = ~paste0(name," (lon:",round(x.centroid,2),", lat:",round(y.centroid,2),")"),
              labelOptions = labelOptions(textOnly = TRUE))

Durata delle soste in porto

Recuperiamo i dati sulle pause delle navi, per vedere quanto si fermano in ciascun porto. Vogliamo associare anche le soste brevi ai “porti” che abbiamo definito prima. Innanzitutto, calcoliamo un raggio di influenza per ciascun porto.

centroids <- dat_cluster %>% 
  group_by(cluster,x.centroid,y.centroid,name) %>% 
  summarize(xsd=sd(x),ysd=sd(y),xm=mean(x),ym=mean(y))
centroids$xsd[is.na(centroids$xsd)]<- 0
centroids$ysd[is.na(centroids$ysd)]<- 0
centroids$radius <- pmax(5, spDists(cbind(centroids$xm-2*centroids$xsd,
                                          centroids$ym-2*centroids$ysd),
                                    cbind(centroids$xm+2*centroids$xsd,
                                          centroids$ym+2*centroids$ysd), 
                                    longlat = T, diagonal=T))
qq <- round(quantile(centroids$radius,c(0.25,0.75)))
centroids$radius <- round(pmin(qq[2], centroids$radius))
centroids$radius <- pmax(qq[1], centroids$radius)

Associamo ogni nave in sosta al porto più vicino, purché la nave sia in sosta per più di un’ora e sia entro il raggio di influenza del porto.

dd <- spDists(as.matrix(dat_pause[,c("x","y")]), 
              as.matrix(centroids[,c("x.centroid","y.centroid")]),
              longlat = T, diagonal=F)
dat_pause$cluster <- centroids$cluster[apply(X = dd,MARGIN = 1,FUN = which.min)]
dat_pause$from_centroid <- apply(X = dd,MARGIN = 1,FUN = min)
dat_pause <- merge(dat_pause,
                   centroids %>% ungroup() %>% select(cluster,radius,name), 
                   by.x="cluster", by.y="cluster")
idx <- dat_pause$from_centroid>dat_pause$radius|dat_pause$duration<60
dat_pause$cluster[idx] <- NA
dat_pause$name[idx] <- NA
leaflet(data = dat_pause %>%
          filter(!is.na(cluster))) %>% 
  addTiles() %>%
  addCircleMarkers(~x, ~y, fillOpacity = 1,
                   color = ~pal(as.numeric(cluster)), radius=3,
                   stroke = FALSE, label = ~name,
                   labelOptions = labelOptions(textOnly = TRUE))
dat_pause %>%
  mutate(dur_h=as.numeric(duration)/60) %>%
  group_by(name) %>% 
  summarize(dur_mean=round(mean(dur_h),1),
            dur_median=round(median(dur_h),1),
            dur_p25=round(quantile(dur_h,0.25),1),
            dur_p75=round(quantile(dur_h,0.75)),
            dur_max=round(max(dur_h)),
            n=n()) %>%
  arrange(desc(n)) -> port_pause
  DT::datatable(port_pause, filter="top")

Percorsi delle navi

Possiamo disegnare i percorsi delle navi. In questo caso conviene selezionare solo alcune giornate.

plot_daily_routes(dat %>% mutate(day=as.Date(substr(date,1,10))) %>%
                    filter(day>as.Date("2018-01-28"),
                           day<as.Date("2018-02-02")))

Oppure, rappresentiamo la densità delle traiettorie, selezionando un periodo più lungo.

traj_levels(dat %>% mutate(day=as.Date(substr(date,1,10))) %>%
                    filter(day>as.Date("2018-01-28"),
                           day<as.Date("2018-02-12")))