# Processing big geographic data sets with R

#### Manuel Campagnolo (ISA/ULisboa)

#### September 12-14, 2018

# 1 Packages

    library(raster)
    library(rgdal)
    library(sf)  # new package for geographic vector datasets (replaces sp)
    library(mapview)  #mapView
    library(RANN)  # nn2
    library(data.table)  # fread
    library(fasterize)
    library(tidyverse)  #  contains dplyr
    library(knitr)
    library(date)
    library(igraph)
    library(randomcoloR)
    library(concaveman)

# 2 Problem: determining meaningful spatial-temporal patches of active fires

Let’s recall the approximate distribution of satellite (VIIRS) active fire observations over Portugal in 2017: it gives some idea of the spatial distribution of fires but doesn’t distinguish dates. The goal of the problem at hand is to determine meaningful spatial-temporal patches (or clusters) of active fires.

    af <- data.table::fread(file.path(getwd(), "datasets", "viirs-fires-2017-world.csv"))
    afpt <- af[latitude > 36 & latitude < 43 & longitude < -6 & longitude > -10]  # Continental Portugal
    latlon <- CRS("+proj=longlat +datum=WGS84")
    y <- raster(crs = latlon, resolution = 0.01, xmn = -10, xmx = -6, ymn = 36, ymx = 42)
    r <- rasterize(x = data.frame(lon = afpt$longitude, lat = afpt$latitude), y = y,
        field = afpt$frp, fun = max)
    mapview(r)

The goal is to define spatial-temporal patches of fire activity. Let’s first convert the dates in the active fire data set into Julian dates. For date::as.date the starting day is January 1, 1960, although this is irrelevant for our example. One could additionally take into account the time of acquisition available in the file and create a time in hours or minutes for each active fires. We are not going to do that to keep the example as simple as possible.

    dt <- date::as.date(x = afpt$acq_date, order = "ymd") %>% as.integer()
    hist(dt, xlabel = "Julian dates (starting January 1, 1960", ylabel = "number of active fires",
        main = "VIIRS active fires: Portugal 2017")

We have three dimensions: one is time, which is measured in days, and the other two are in space. To measure distances on the surface, one should avoid geographic coordinates. We can re-project the data onto a local coordinate reference system as the official CRS for Portugal (PT-TM06/ETRS89, EPSG:3763).

To make the sequence of operations easier to follow, we use magrittr’s pipe operator for chaining commands.

    sfpt <- as.data.frame(afpt) %>% st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
        st_transform(crs = 3763)

# 3 Determining spatial distances between active fires

Spatial-temporal patches depend on the distance between active fires. Since coordinates are now in meters, we set the parameter d with a value in meters. For the time dimension, we define a time gap t between active fires.

    d <- 2000  # meters
    t <- 2  # days
    K <- 20  # neighbors to consider

Next we build a graph which vertices represent active fires. In the graph, active fires are going to be connected if the distance between them is shorter than d and if the tile gap is not larger than t. We start by determining the list of edges that connect vertices. This is going to be ultimately represented by a data.table with columns v1 and v2, where each pair (v1,v2) represents an edge in the graph.

    knn <- nn2(data = st_coordinates(sfpt), k = K, searchtype = "radius", radius = d)
    head(knn$nn.idx[, 1:10], 5)

         [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
    [1,]    1    2    3 1059    0    0    0    0    0     0
    [2,]    2    1    3 1059    0    0    0    0    0     0
    [3,]    3    2    1 1059    0    0    0    0    0     0
    [4,]    4  669    0    0    0    0    0    0    0     0
    [5,]    5   36   34    8    6   37   35   22    9    21

    head(knn$nn.dists[, 1:5], 3)

         [,1]     [,2]     [,3]      [,4]          [,5]
    [1,]    0 690.2998 938.2042 1081.2232 1.340781e+154
    [2,]    0 690.2998 723.8893  977.0478 1.340781e+154
    [3,]    0 723.8893 938.2042 1680.1948 1.340781e+154

Note that the distance to the nearest neighbor is always 0. This is because in general the nearest neighbor is the vertex itself (if there are two vertices with precisely the same coordinates, the nearest neighbor can be a distinct vertex, but the distance is still 0). The next step is to determine all the edges between neighbor observations. In order to do this we rearrange the columns of knn$nn.idx. The first command is to create a first column with the indices of all observations (1,2,3…). The next command moves down to the second column of edges all neighbors, so the resulting matrix edges has just two columns, where the first column has the indices of the observations (repeated K=ncol(neigh) times), and the second column has the respective j-th neighbor.

    neigh <- cbind(1:nrow(knn$nn.idx), knn$nn.idx)
    head(neigh[, 1:10], 4)

         [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
    [1,]    1    1    2    3 1059    0    0    0    0     0
    [2,]    2    2    1    3 1059    0    0    0    0     0
    [3,]    3    3    2    1 1059    0    0    0    0     0
    [4,]    4    4  669    0    0    0    0    0    0     0

    edges <- data.table(v1 = rep(neigh[, 1], ncol(neigh) - 1), v2 = as.vector(neigh[,
        -1]))
    dim(edges)

    [1] 519180      2

At this point, one should remove pairs where one vertex is missing (in that case the index in the second column is 0).

    edges[v2 == 0]  # non edges

              v1 v2
        1:    28  0
        2:    29  0
        3:    60  0
        4:    64  0
        5:    92  0
       ---
    97198: 25943  0
    97199: 25948  0
    97200: 25949  0
    97201: 25953  0
    97202: 25959  0

    edges <- edges[v2 != 0]  # remove non edges

# 4 Adding temporal constraints

So far, we only used spatial coordinates, but we want our patches to satisfy a temporal restriction as well. As mentioned above, two active fires should be connected only if their acquisition dates differs at most by t days. In order to enforce this, we first need to add to our matrix the acquisition dates available in vector dt. Note, and this is a key aspect of the construction, that observations are indexed as in the original data set (the same order and length). Therefore the index of each observation in edges is precisely the same index than in dt.

    edges[, `:=`(d1, dt[v1])]  # acquisition dates
    edges[, `:=`(d2, dt[v2])]  # acquisition dates
    edges <- edges[abs(d1 - d2) <= t]

Just for sake of clarity, we also add the information about the coordinates of the observations, although this is not necessary for the task at hand.

    xy <- st_coordinates(sfpt)  # for clarity, but optional
    edges[, `:=`(x1, xy[v1, "X"])]  # coordinates, optional
    edges[, `:=`(y1, xy[v1, "Y"])]  # coordinates, optional
    edges[, `:=`(x2, xy[v2, "X"])]  # coordinates, optional
    edges[, `:=`(y2, xy[v2, "Y"])]  # coordinates, optional
    edges[1:3, ]

       v1 v2    d1    d2       x1       y1       x2       y2
    1:  1  1 20820 20820 95620.27 331186.0 95620.27 331186.0
    2:  2  2 20820 20820 95569.88 330497.5 95569.88 330497.5
    3:  3  3 20820 20820 96292.95 330532.0 96292.95 330532.0

Let’s double check that all distances are within d and that all date gaps are within t:

    max(edges[, sqrt((x1 - x2)^2 + (y1 - y2)^2)])

    [1] 1999.929

    max(edges[, abs(d1 - d2)])

    [1] 2

    sum(duplicated(edges))  # check that there are no duplicated edges

    [1] 0

# 5 Building and exploring the graph

We have now a list of 354232 and we can build a graph that represents active fires and their proximity relations. There are many types of graphs, including directed and undirected graphs. In our case, we only need an undirected graph, where the directions on the edges are not relevant. Note that the construction above doesn’t guarantee that there are no multiple edges (i.e. two edges (v1,v2) and (v2,v1)). We could remove those edges directly from data.table edges but it is better to use available functions from package igraph which are properly designed for that goal. In our case, we use function igraph::simplify that removes loops and multiple edges.

An introduction about drawing graphs with igraph can be found here.

    Gall <- igraph::graph_from_data_frame(d = edges, directed = FALSE)
    G <- igraph::simplify(Gall)  # removes loops and multiple edges

Let’s just confirm that the names of the vertices in the graph correspond to the indices of the vertices in the table edges.

    v <- as.integer(vertex_attr(G)$name)  # vertex indices which match row indices of edges
    range(v) == range(edges[, .(v1, v2)])  # confirm index numering

    [1] TRUE TRUE

Next, we use a igraph functions that returns the connected components in G. Those are the clusters of vertices that are connected in the graph. Notice that the result is a named vector with one value per vertex: the name is the name of the vertex (of character type) and the value is the component number. This tell us what is the component for each vertex.

    CL <- igraph::components(G)
    head(CL$membership)  # returns vertice name and vertice cluster

    1 2 3 4 5 6
    1 1 1 2 3 3

    head(names(CL$membership))  # names (character)

    [1] "1" "2" "3" "4" "5" "6"

Let’s consider one component of the graph and see what it looks like. It contains all vertices in that component and the edges between them. Function induced_graph returns the sub-graph with a subset of edges and all connections between them.

    head(CL$csize)

    [1]  3  1 14  2  5  1

    idxcomp <- 2500
    # vertices of the component
    vids <- as.integer(names(CL$membership)[CL$membership == idxcomp])
    # subgraph
    S <- igraph::induced_subgraph(G, vids = which(CL$membership == idxcomp))

We can plot the graph, and add to the plot additional information using labels and colors. Here, we label each vertex with the active fire acquisition date, and we color each vertex according to the FRP (Fire Radiative Power, expressed in Mega Watts). A vertex attribute (e.g. color) can be set for the graph with set.vertex.attribute. In alternative, the attribute can be set as a plot option.

    # colors for subgraph: depending on FRP
    cores <- colorRampPalette(c(rgb(1, 1, 0), rgb(1, 0, 0)))(length(V(S)))[order(sfpt$frp[vids])]
    S <- set.vertex.attribute(S, "color", value = cores)
    # color is a new attribute
    vertex.attributes(S)

    $name
     [1] "10777" "10778" "11165" "11166" "11167" "11168" "11169" "11170" "11171"
    [10] "11172" "11173" "11175" "11176" "11177" "11611"

    $color
     [1] "#FFDA00" "#FFB600" "#FF0000" "#FF5B00" "#FF4800" "#FF7F00" "#FF6D00"
     [8] "#FFC800" "#FFA300" "#FF9100" "#FF3600" "#FF1200" "#FF2400" "#FFFF00"
    [15] "#FFEC00"

    # plot using xy coordinates, create layout for ploting
    lxy <- st_coordinates(sfpt)[vids, ]

Plotting defining labels and colors for the vertices, and lxy layout for the vertices.

    plot(S, layout = lxy, vertex.label = sfpt$acq_date[vids], vertex.size = 30, vertex.color = cores)
    plot(S, layout = lxy, vertex.label = sfpt$acq_time[vids], vertex.size = 30, vertex.color = cores)

# 6 Mapping spatial-temporal clusters of active fires

Let’s go back to our sf object and add the information about the clusters in sfpt. This is very easy to do since it is just a new column with the cluster membership. Note that this works because the graph G was originally derived from the table edges, which, in turn, was build from all points in sfpt. Therefore, the vertex’ indices match the rows of sfpt. For sake of completude, we re-compute components from graph G and then we incorporate that information in the spatial object sfpt.

    CL <- igraph::components(G)
    sfpt$clust <- CL$membership

Let’s then make the first map of spatial temporal patches. Here we assign the same color to all points that belong to the same patch, so we can visualize them in the map. We use random colors to make it unlikely that two neighbor patches have similar colors.

    n <- 20
    palette <- distinctColorPalette(n)
    cores <- rep(palette, times = max(sfpt$clust)/n)
    mapview(sfpt, color = cores[sfpt$clust], cex = 4)