Motivation

https://oscarperpinan.github.io/rastervis/

https://oscarperpinan.github.io/rastervis/

  • Aim: recreate a “vectorplot” in ggplot2, but using a real dataset rather than a synthetic example
  • Specifically -
    • We have a dataset with storm trajectories
    • We want to display a regular grid with of arrow segments, with direction and length representing the direction and intensity of storms passing through that grid cell
  • How do we do that?
    • Convert the trajectories to individual segments
    • Calculate the azimuth of each storm segment
    • Make a regular grid
    • Summarize the average azimuth and total length of storm segments
    • Calculate gridded segments, according average azimuth and total length in each cell
    • Plot

Storms dataset

library(dplyr)

storms = as.data.frame(storms)
head(storms)
##   name year month day hour  lat  long              status category wind
## 1  Amy 1975     6  27    0 27.5 -79.0 tropical depression       -1   25
## 2  Amy 1975     6  27    6 28.5 -79.0 tropical depression       -1   25
## 3  Amy 1975     6  27   12 29.5 -79.0 tropical depression       -1   25
## 4  Amy 1975     6  27   18 30.5 -79.0 tropical depression       -1   25
## 5  Amy 1975     6  28    0 31.5 -78.8 tropical depression       -1   25
## 6  Amy 1975     6  28    6 32.4 -78.7 tropical depression       -1   25
##   pressure ts_diameter hu_diameter
## 1     1013          NA          NA
## 2     1013          NA          NA
## 3     1013          NA          NA
## 4     1013          NA          NA
## 5     1012          NA          NA
## 6     1012          NA          NA

Storm ID

  • There are storms with the same name in different years
  • To identify individual storm tracks we add a name_year ID variable -
storms$name_year = paste(storms$name, storms$year)
head(storms)
##   name year month day hour  lat  long              status category wind
## 1  Amy 1975     6  27    0 27.5 -79.0 tropical depression       -1   25
## 2  Amy 1975     6  27    6 28.5 -79.0 tropical depression       -1   25
## 3  Amy 1975     6  27   12 29.5 -79.0 tropical depression       -1   25
## 4  Amy 1975     6  27   18 30.5 -79.0 tropical depression       -1   25
## 5  Amy 1975     6  28    0 31.5 -78.8 tropical depression       -1   25
## 6  Amy 1975     6  28    6 32.4 -78.7 tropical depression       -1   25
##   pressure ts_diameter hu_diameter name_year
## 1     1013          NA          NA  Amy 1975
## 2     1013          NA          NA  Amy 1975
## 3     1013          NA          NA  Amy 1975
## 4     1013          NA          NA  Amy 1975
## 5     1012          NA          NA  Amy 1975
## 6     1012          NA          NA  Amy 1975

To points

  • The table is converted to a point layer with st_as_sf using the long and lat columns -
library(sf)

pnt = st_as_sf(storms, coords = c("long", "lat"), crs = 4326)
pnt
## Simple feature collection with 10010 features and 12 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -6 ymax: 51.9
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##    name year month day hour              status category wind pressure
## 1   Amy 1975     6  27    0 tropical depression       -1   25     1013
## 2   Amy 1975     6  27    6 tropical depression       -1   25     1013
## 3   Amy 1975     6  27   12 tropical depression       -1   25     1013
## 4   Amy 1975     6  27   18 tropical depression       -1   25     1013
## 5   Amy 1975     6  28    0 tropical depression       -1   25     1012
## 6   Amy 1975     6  28    6 tropical depression       -1   25     1012
## 7   Amy 1975     6  28   12 tropical depression       -1   25     1011
## 8   Amy 1975     6  28   18 tropical depression       -1   30     1006
## 9   Amy 1975     6  29    0      tropical storm        0   35     1004
## 10  Amy 1975     6  29    6      tropical storm        0   40     1002
##    ts_diameter hu_diameter name_year           geometry
## 1           NA          NA  Amy 1975   POINT (-79 27.5)
## 2           NA          NA  Amy 1975   POINT (-79 28.5)
## 3           NA          NA  Amy 1975   POINT (-79 29.5)
## 4           NA          NA  Amy 1975   POINT (-79 30.5)
## 5           NA          NA  Amy 1975 POINT (-78.8 31.5)
## 6           NA          NA  Amy 1975 POINT (-78.7 32.4)
## 7           NA          NA  Amy 1975   POINT (-78 33.3)
## 8           NA          NA  Amy 1975     POINT (-77 34)
## 9           NA          NA  Amy 1975 POINT (-75.8 34.4)
## 10          NA          NA  Amy 1975   POINT (-74.8 34)

Plot

plot(pnt, max.plot = 12)

Points to lines

  • Next we want to create a layer of line segments between all consecutive points of each storm
  • First we split the point layer by storm ID, to get a list of POINT subsets -
lines = split(pnt, pnt$name_year)
lines[1]
## $`AL011993 1993`
## Simple feature collection with 8 features and 12 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -84 ymin: 21.5 xmax: -71.8 ymax: 27.8
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##          name year month day hour              status category wind
## 2990 AL011993 1993     5  31   12 tropical depression       -1   25
## 2991 AL011993 1993     5  31   18 tropical depression       -1   25
## 2992 AL011993 1993     6   1    0 tropical depression       -1   25
## 2993 AL011993 1993     6   1    6 tropical depression       -1   25
## 2994 AL011993 1993     6   1   12 tropical depression       -1   30
## 2995 AL011993 1993     6   1   18 tropical depression       -1   30
## 2996 AL011993 1993     6   2    0 tropical depression       -1   30
## 2997 AL011993 1993     6   2    6 tropical depression       -1   30
##      pressure ts_diameter hu_diameter     name_year           geometry
## 2990     1003          NA          NA AL011993 1993   POINT (-84 21.5)
## 2991     1002          NA          NA AL011993 1993   POINT (-82 22.3)
## 2992     1000          NA          NA AL011993 1993 POINT (-80.3 23.2)
## 2993     1000          NA          NA AL011993 1993   POINT (-79 24.5)
## 2994      999          NA          NA AL011993 1993 POINT (-77.5 25.4)
## 2995      999          NA          NA AL011993 1993 POINT (-75.8 26.1)
## 2996      999          NA          NA AL011993 1993   POINT (-74 26.7)
## 2997      999          NA          NA AL011993 1993 POINT (-71.8 27.8)
  • Then, we combine the points to MULTIPOINT with st_combine -
lines = lapply(lines, st_combine)
lines[1]
## $`AL011993 1993`
## Geometry set for 1 feature 
## geometry type:  MULTIPOINT
## dimension:      XY
## bbox:           xmin: -84 ymin: 21.5 xmax: -71.8 ymax: 27.8
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
  • Casting to LINESTRING -
lines = lapply(lines, st_cast, "LINESTRING")
lines[1]
## $`AL011993 1993`
## Geometry set for 1 feature 
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -84 ymin: 21.5 xmax: -71.8 ymax: 27.8
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
  • The list of LINESTRING can be combined back to a geometry column with c -
geometry = do.call(c, lines)
geometry
## Geometry set for 426 features 
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -6 ymax: 51.9
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 5 geometries:
  • The corresponding attribute name_year is kept in the list element names -
attr = data.frame(name_year = names(lines), stringsAsFactors = FALSE)
head(attr)
##       name_year
## 1 AL011993 1993
## 2 AL012000 2000
## 3 AL021992 1992
## 4 AL021994 1994
## 5 AL021999 1999
## 6 AL022000 2000
  • Combining the geometry and attributes back to an sf layer -
lines = st_sf(geometry, attr)
lines
## Simple feature collection with 426 features and 1 field
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -6 ymax: 51.9
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##        name_year                       geometry
## 1  AL011993 1993 LINESTRING (-84 21.5, -82 2...
## 2  AL012000 2000 LINESTRING (-93 21, -92.8 2...
## 3  AL021992 1992 LINESTRING (-85.5 24.5, -85...
## 4  AL021994 1994 LINESTRING (-78.9 32.2, -79...
## 5  AL021999 1999 LINESTRING (-95 20.2, -96.3...
## 6  AL022000 2000 LINESTRING (-19.8 9.5, -21 ...
## 7  AL022001 2001 LINESTRING (-42.1 10.9, -43...
## 8  AL022003 2003 LINESTRING (-40.8 9.5, -42....
## 9  AL022006 2006 LINESTRING (-66.4 39.1, -65...
## 10 AL031987 1987 LINESTRING (-93.6 26.3, -94...

Lines

plot(lines)

Lines to segments

  • To calculate segment azimuth, we split the lines to segments -
library(shadow)

seg = 
  lines %>% 
  as("Spatial") %>% 
  toSeg %>% 
  st_as_sf
seg
## Simple feature collection with 9584 features and 1 field
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -6 ymax: 51.9
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##         name_year                       geometry
## 1 1 AL011993 1993 LINESTRING (-84 21.5, -82 2...
## 1 2 AL011993 1993 LINESTRING (-82 22.3, -80.3...
## 1 3 AL011993 1993 LINESTRING (-80.3 23.2, -79...
## 1 4 AL011993 1993 LINESTRING (-79 24.5, -77.5...
## 1 5 AL011993 1993 LINESTRING (-77.5 25.4, -75...
## 1 6 AL011993 1993 LINESTRING (-75.8 26.1, -74...
## 1 7 AL011993 1993 LINESTRING (-74 26.7, -71.8...
## 2 1 AL012000 2000 LINESTRING (-93 21, -92.8 2...
## 2 2 AL012000 2000 LINESTRING (-92.8 20.9, -93...
## 2 3 AL012000 2000 LINESTRING (-93.1 20.7, -93...
  • Casting to MULTIPOINT; each segment is converted to a pair of points -
seg_pnt = st_cast(seg, "MULTIPOINT")
seg_pnt
## Simple feature collection with 9584 features and 1 field
## geometry type:  MULTIPOINT
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -6 ymax: 51.9
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##        name_year                       geometry
## 1  AL011993 1993 MULTIPOINT (-84 21.5, -82 2...
## 2  AL011993 1993 MULTIPOINT (-82 22.3, -80.3...
## 3  AL011993 1993 MULTIPOINT (-80.3 23.2, -79...
## 4  AL011993 1993 MULTIPOINT (-79 24.5, -77.5...
## 5  AL011993 1993 MULTIPOINT (-77.5 25.4, -75...
## 6  AL011993 1993 MULTIPOINT (-75.8 26.1, -74...
## 7  AL011993 1993 MULTIPOINT (-74 26.7, -71.8...
## 8  AL012000 2000 MULTIPOINT (-93 21, -92.8 2...
## 9  AL012000 2000 MULTIPOINT (-92.8 20.9, -93...
## 10 AL012000 2000 MULTIPOINT (-93.1 20.7, -93...

Segment azimuth

  • Azimuth can be calculated with bearingRhumb from package geosphere
  • For example, here are the coordinates of the first storm segment -
x = st_coordinates(seg_pnt[1, ])
x
##        X    Y L1
## [1,] -84 21.5  1
## [2,] -82 22.3  1
  • And here is its azimuth -
library(geosphere)

bearingRhumb(x[1, 1:2], x[2, 1:2])
## [1] 66.67839
  • Now, the same can be done for all segments
  • We create a list of segment coordinates -
coords = lapply(st_geometry(seg_pnt), st_coordinates)
head(coords)
## [[1]]
##        X    Y L1
## [1,] -84 21.5  1
## [2,] -82 22.3  1
## 
## [[2]]
##          X    Y L1
## [1,] -82.0 22.3  1
## [2,] -80.3 23.2  1
## 
## [[3]]
##          X    Y L1
## [1,] -80.3 23.2  1
## [2,] -79.0 24.5  1
## 
## [[4]]
##          X    Y L1
## [1,] -79.0 24.5  1
## [2,] -77.5 25.4  1
## 
## [[5]]
##          X    Y L1
## [1,] -77.5 25.4  1
## [2,] -75.8 26.1  1
## 
## [[6]]
##          X    Y L1
## [1,] -75.8 26.1  1
## [2,] -74.0 26.7  1
  • Then apply bearingRhumb on the list -
seg$az = sapply(coords, function(x) bearingRhumb(x[1, 1:2], x[2, 1:2]))
seg = seg[!is.na(seg$az), ] # Remove zero distance segments
seg
## Simple feature collection with 9548 features and 2 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -6 ymax: 51.9
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##         name_year                       geometry        az
## 1 1 AL011993 1993 LINESTRING (-84 21.5, -82 2...  66.67839
## 1 2 AL011993 1993 LINESTRING (-82 22.3, -80.3...  60.14064
## 1 3 AL011993 1993 LINESTRING (-80.3 23.2, -79...  42.44541
## 1 4 AL011993 1993 LINESTRING (-79 24.5, -77.5...  56.50471
## 1 5 AL011993 1993 LINESTRING (-77.5 25.4, -75...  65.43174
## 1 6 AL011993 1993 LINESTRING (-75.8 26.1, -74...  69.58745
## 1 7 AL011993 1993 LINESTRING (-74 26.7, -71.8...  60.64520
## 2 1 AL012000 2000 LINESTRING (-93 21, -92.8 2... 118.16430
## 2 2 AL012000 2000 LINESTRING (-92.8 20.9, -93... 234.50559
## 2 3 AL012000 2000 LINESTRING (-93.1 20.7, -93... 284.96749
  • We also need the segment centroids for later -
seg_ctr = st_centroid(seg)
seg_ctr
## Simple feature collection with 9548 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -109.2 ymin: 7.3 xmax: -6.55 ymax: 51.65
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##         name_year             geometry        az
## 1 1 AL011993 1993     POINT (-83 21.9)  66.67839
## 1 2 AL011993 1993 POINT (-81.15 22.75)  60.14064
## 1 3 AL011993 1993 POINT (-79.65 23.85)  42.44541
## 1 4 AL011993 1993 POINT (-78.25 24.95)  56.50471
## 1 5 AL011993 1993 POINT (-76.65 25.75)  65.43174
## 1 6 AL011993 1993   POINT (-74.9 26.4)  69.58745
## 1 7 AL011993 1993  POINT (-72.9 27.25)  60.64520
## 2 1 AL012000 2000  POINT (-92.9 20.95) 118.16430
## 2 2 AL012000 2000  POINT (-92.95 20.8) 234.50559
## 2 3 AL012000 2000  POINT (-93.3 20.75) 284.96749
plot(seg_ctr)

Make grid

  • Making a regular grid with a unique ID -
grid = st_make_grid(lines, cellsize = 2)
grid = st_sf(grid, data.frame(id = 1:length(grid)))
grid
## Simple feature collection with 1196 features and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -5.3 ymax: 53.2
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##    id                           grid
## 1   1 POLYGON ((-109.3 7.2, -107....
## 2   2 POLYGON ((-107.3 7.2, -105....
## 3   3 POLYGON ((-105.3 7.2, -103....
## 4   4 POLYGON ((-103.3 7.2, -101....
## 5   5 POLYGON ((-101.3 7.2, -99.3...
## 6   6 POLYGON ((-99.3 7.2, -97.3 ...
## 7   7 POLYGON ((-97.3 7.2, -95.3 ...
## 8   8 POLYGON ((-95.3 7.2, -93.3 ...
## 9   9 POLYGON ((-93.3 7.2, -91.3 ...
## 10 10 POLYGON ((-91.3 7.2, -89.3 ...
plot(grid)

Trajectory length per grid cell

  • Next we need to find the total trajectory length per grid cell
  • First we intersect the grid and seg layers, and calculate the length per segment -
grid_seg = st_intersection(grid, seg)
grid_seg$length = st_length(grid_seg)
grid_seg
## Simple feature collection with 17574 features and 4 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -6 ymax: 51.9
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##        id     name_year       az                           grid
## 377   377 AL011993 1993 66.67839 LINESTRING (-84 21.5, -83.3...
## 378   378 AL011993 1993 66.67839 LINESTRING (-83.3 21.78, -8...
## 378.1 378 AL011993 1993 60.14064 LINESTRING (-82 22.3, -81.3...
## 379   379 AL011993 1993 60.14064 LINESTRING (-81.3 22.67059,...
## 431   431 AL011993 1993 60.14064             POINT (-80.3 23.2)
## 379.1 379 AL011993 1993 42.44541             POINT (-80.3 23.2)
## 431.1 431 AL011993 1993 42.44541 LINESTRING (-80.3 23.2, -79...
## 432   432 AL011993 1993 42.44541 LINESTRING (-79.3 24.2, -79...
## 432.1 432 AL011993 1993 56.50471 LINESTRING (-79 24.5, -77.8...
## 484   484 AL011993 1993 56.50471 LINESTRING (-77.83333 25.2,...
##            length
## 377    78817.99 m
## 378   146033.02 m
## 378.1  82904.18 m
## 379   118144.19 m
## 431        0.00 m
## 379.1      0.00 m
## 431.1 150556.21 m
## 432    45065.53 m
## 432.1 141126.12 m
## 484    40220.19 m
  • Then we aggregate the table, summarizing total length -
grid_seg_agg = aggregate(
  x = st_set_geometry(grid_seg[, "length"], NULL), 
  by = data.frame(id = grid_seg$id), 
  FUN = sum
)
head(grid_seg_agg)
##   id     length
## 1 28  49221.57 
## 2 30 109896.87 
## 3 31 159778.13 
## 4 32 198491.48 
## 5 33 108670.34 
## 6 34 238682.00
  • And join back to the grid -
grid = left_join(grid, grid_seg_agg, "id")
grid
## Simple feature collection with 1196 features and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -5.3 ymax: 53.2
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##    id length                           grid
## 1   1    NA  POLYGON ((-109.3 7.2, -107....
## 2   2    NA  POLYGON ((-107.3 7.2, -105....
## 3   3    NA  POLYGON ((-105.3 7.2, -103....
## 4   4    NA  POLYGON ((-103.3 7.2, -101....
## 5   5    NA  POLYGON ((-101.3 7.2, -99.3...
## 6   6    NA  POLYGON ((-99.3 7.2, -97.3 ...
## 7   7    NA  POLYGON ((-97.3 7.2, -95.3 ...
## 8   8    NA  POLYGON ((-95.3 7.2, -93.3 ...
## 9   9    NA  POLYGON ((-93.3 7.2, -91.3 ...
## 10 10    NA  POLYGON ((-91.3 7.2, -89.3 ...
plot(grid)

Average segment azimuth

  • Azimuth should be averaged using circular statistics, for example -
library(circular)

angles = c(2, 359)

# Simple mean
mean(angles)
## [1] 180.5
# Circular mean
anglecir = circular(angles, type="angles", units="degrees", modulo="2pi", template="geographics")
as.numeric(mean(anglecir))
## [1] 0.5
  • And now for our data, averaging individual segment azimuths per grid cell -
f = function(x) {
  anglecir = circular(x, type="angles", units="degrees", modulo="2pi", template="geographics")
  m = mean(anglecir)
  as.numeric(m)
}
grid_b = aggregate(x = seg_ctr[, "az"], by = grid, FUN = f)
grid_b
## Simple feature collection with 1196 features and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -109.3 ymin: 7.2 xmax: -5.3 ymax: 53.2
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##    az                       geometry
## 1  NA POLYGON ((-109.3 7.2, -107....
## 2  NA POLYGON ((-107.3 7.2, -105....
## 3  NA POLYGON ((-105.3 7.2, -103....
## 4  NA POLYGON ((-103.3 7.2, -101....
## 5  NA POLYGON ((-101.3 7.2, -99.3...
## 6  NA POLYGON ((-99.3 7.2, -97.3 ...
## 7  NA POLYGON ((-97.3 7.2, -95.3 ...
## 8  NA POLYGON ((-95.3 7.2, -93.3 ...
## 9  NA POLYGON ((-93.3 7.2, -91.3 ...
## 10 NA POLYGON ((-91.3 7.2, -89.3 ...
plot(grid_b)

Grid centroids

  • For the final visualization we need the grid cell centroids -
grid_ctr = st_centroid(grid)
plot(grid_ctr)

  • The average azimuths can be joined to the grid centroids using spatial join with st_join -
grid_ctr = st_join(grid_ctr, grid_b)
grid_ctr = grid_ctr[!is.na(grid_ctr$az) & !is.na(grid_ctr$length), ]
grid_ctr
## Simple feature collection with 690 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -108.3 ymin: 8.2 xmax: -6.3 ymax: 52.2
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##    id    length       az              grid
## 30 30 109896.9  270.0000 POINT (-50.3 8.2)
## 31 31 159778.1  276.7965 POINT (-48.3 8.2)
## 32 32 198491.5  273.8620 POINT (-46.3 8.2)
## 33 33 108670.3  286.1382 POINT (-44.3 8.2)
## 34 34 238682.0  283.9600 POINT (-42.3 8.2)
## 35 35 178224.4  281.4519 POINT (-40.3 8.2)
## 36 36 220176.8  271.7039 POINT (-38.3 8.2)
## 37 37 221128.4  262.3118 POINT (-36.3 8.2)
## 38 38 293552.9  276.6851 POINT (-34.3 8.2)
## 39 39 213754.8  274.1292 POINT (-32.3 8.2)
plot(grid_ctr)

Shift grid centroids

  • Grid cell centroids will comprise the segment staring point
  • To get the ending point, we shift the starting point -
    • In the direction of average storm segment azimuth
    • With distance proportional to total segment length
library(scales)

grid_ctr_shifted = 
  grid_ctr %>% 
  as("Spatial") %>% 
  shiftAz(
    az = grid_ctr$az, 
    dist = rescale(as.numeric(grid_ctr$length), c(0, 3))
  ) %>% 
  st_as_sf
plot(st_geometry(grid_ctr[1:3, ]))
plot(st_geometry(grid_ctr_shifted[1:3, ]), col = "red", add = TRUE)

Final table

  • The final table -
dat = cbind(
  data.frame(
    lon0 = st_coordinates(grid_ctr)[, 1],
    lat0 = st_coordinates(grid_ctr)[, 2],
    length = grid_ctr$length,
    az = grid_ctr$az
  ),
  data.frame(
    lon1 = st_coordinates(grid_ctr_shifted)[, 1],
    lat1 = st_coordinates(grid_ctr_shifted)[, 2]
  )
)
dat = dat[, c("lon0", "lat0", "lon1", "lat1", "length", "az")]
head(dat)
##    lon0 lat0      lon1     lat1   length       az
## 1 -50.3  8.2 -50.33352 8.200000 109896.9 270.0000
## 2 -48.3  8.2 -48.35888 8.207017 159778.1 276.7965
## 3 -46.3  8.2 -46.37912 8.205341 198491.5 273.8620
## 4 -44.3  8.2 -44.33159 8.209141 108670.3 286.1382
## 5 -42.3  8.2 -42.39711 8.224139 238682.0 283.9600
## 6 -40.3  8.2 -40.36745 8.213665 178224.4 281.4519

Initial plot

library(ggplot2)

ggplot(dat) +
  geom_sf(data = lines, colour = "lightgrey") +
  geom_segment(aes(x = lon0, y = lat0, xend = lon1, yend = lat1)) +
  theme_bw()

Circular colors

  • Setting colors -
cols = circular.colors(n = 12)
cols
##  [1] "#FF0000" "#FF8B00" "#E8FF00" "#5DFF00" "#00FF2E" "#00FFB9" "#00B9FF"
##  [8] "#002EFF" "#5D00FF" "#E800FF" "#FF008B" "#FF0000"
  • Setting color breaks -
dat$b1 = cut(dat$az, breaks = seq(0, 360, 30))
head(dat)
##    lon0 lat0      lon1     lat1   length       az        b1
## 1 -50.3  8.2 -50.33352 8.200000 109896.9 270.0000 (240,270]
## 2 -48.3  8.2 -48.35888 8.207017 159778.1 276.7965 (270,300]
## 3 -46.3  8.2 -46.37912 8.205341 198491.5 273.8620 (270,300]
## 4 -44.3  8.2 -44.33159 8.209141 108670.3 286.1382 (270,300]
## 5 -42.3  8.2 -42.39711 8.224139 238682.0 283.9600 (270,300]
## 6 -40.3  8.2 -40.36745 8.213665 178224.4 281.4519 (270,300]

Final plot

ggplot(dat) +
  geom_sf(data = lines, colour = "lightgrey", alpha = 0.1) +
  geom_segment(
    aes(x = lon0, y = lat0, xend = lon1, yend = lat1, size = length, colour = b1), 
    arrow = arrow(length = unit(0.015, "npc"))
  ) +
  scale_size_continuous(range = c(0.15, 1), guide = FALSE) +
  scale_color_manual("Azimuth (°)", values = cols) +
  theme_bw() +
  theme(
    axis.title = element_blank()
  )

Summary

  • Often, much of the time in visualization is spent on reshaping the dataset to the right form

Other ideas

  • Equal-area grid
  • Weight by wind speed too, not just length
  • Segment centred on grid cell

https://earth.nullschool.net/