{"id":547,"date":"2018-06-14T20:14:15","date_gmt":"2018-06-15T00:14:15","guid":{"rendered":"https:\/\/micah.waldste.in\/?p=547"},"modified":"2018-06-24T20:15:03","modified_gmt":"2018-06-25T00:15:03","slug":"heatmaps-of-spherical-densities-in-r","status":"publish","type":"post","link":"https:\/\/micah.waldste.in\/blog\/blog\/2018\/06\/heatmaps-of-spherical-densities-in-r\/","title":{"rendered":"Heatmaps of Spherical Densities in R"},"content":{"rendered":"

DISCLAIMER: While I know a thing or two, there's a reasonable chance I got\r\nsome things wrong or at very least there are certainly more efficient ways to\r\ngo about things. Feedback always appreciated!<\/em><\/p>\r\n\r\n

Last time<\/a>\r\nwe made contour maps of densities of points on a globe, now it is\r\ntime to take another step and make heatmaps. We created all the data we needed\r\nwhen creating the contours, but heatmaps add new challenges of dealing with\r\nlarge amounts of raster and polygon data. Lets get to it.<\/p>\r\n\r\n

Set-Up<\/h1>\r\n\r\n

First, we'll make use of a number of libraries and setup our plotting\r\nenvironment:<\/p>\r\n\r\n

library(rgdal)       # For coordinate transforms\r\nlibrary(sp)          # For plotting grid images\r\nlibrary(sf)\r\nlibrary(lwgeom)\r\nlibrary(Directional) # For spherical density functions\r\nlibrary(spData)      # worldmap\r\nlibrary(raster)\r\n<\/code><\/pre>\r\n\r\n
We'll also use the same vmf_density_grid<\/code> function we introduced in the Intro\r\npost.<\/p>\r\n\r\n
vmf_density_grid <- function(u, ngrid = 100) {\r\n  # Translate to (0,180) and (0,360)\r\n  u[,1] <- u[,1] + 90\r\n  u[,2] <- u[,2] + 180\r\n  res <- vmf.kerncontour.new(u, thumb = "none", ret.all = T, full = T,\r\n                             ngrid = ngrid)\r\n\r\n  # Translate back to (-90, 90) and (-180, 180) and create a grid of\r\n  # coordinates\r\n  ret <- expand.grid(Lat = res$Lat - 90, Long = res$Long - 180)\r\n  ret$Density = c(res$d)\r\n  ret\r\n}\r\n<\/code><\/pre>\r\n\r\n
Global Earthquakes Again<\/h2>\r\n\r\n
Global Earthquakes from Northern California Earthquake Data\r\nCenter<\/a> is a great dataset we'll\r\ncontinue to use, so we start with a set of quakes since Jan 1, 1950 of\r\nmagnitude 5.9 or higher.<\/p>\r\n\r\n
For all our heatmaps, we'll start the same as we did for contours, calculating\r\nthe density map:<\/p>\r\n\r\n
grid.size = 100\r\nearthquakes <- read.csv(file.path("..", "data", "earthquakes.csv"))\r\nearthquake.densities <- vmf_density_grid(earthquakes[,c("Latitude",\r\n                                                        "Longitude")],\r\n                                         ngrid = grid.size)\r\n<\/code><\/pre>\r\n\r\n
Once we have the densities, we need to coerce them into a spatial format – in\r\nthis case we'll create a SpatialGridDataFrame<\/code>, matching the grid of densities\r\nwe calculated with vmf_density_grid<\/code>.<\/p>\r\n\r\n
density_matrix <- matrix(earthquake.densities$Density, nrow = grid.size)\r\ndensity_matrix <- t(apply(density_matrix, 2, rev))\r\ngridVals <- data.frame(att=as.vector(density_matrix))\r\ngt <- GridTopology(cellcentre.offset = c(-180 + 180 \/ grid.size,\r\n                                         -90 + 90 \/ grid.size),\r\n                   cellsize = c( 360 \/ grid.size, 180 \/ grid.size),\r\n                   cells.dim = c(grid.size, grid.size))\r\nsGDF <- SpatialGridDataFrame(gt,\r\n                             data = gridVals,\r\n                             proj = "+proj=longlat +datum=WGS84 +no_defs")\r\n\r\nplot(sGDF)\r\nplot(gridlines(sGDF), add = TRUE, col = "grey30", alpha = .1)\r\nplot(st_geometry(world), add = TRUE, col = NA, border = "grey")\r\n<\/code><\/pre>\r\n\r\n
\"plot<\/p>\r\n\r\n
Great, we have a heatmap! But it is in rectangular coordinates, we want to\r\nproject it to something nicer, like a Winkel triple. There's a problem\r\nthough… We can't just re-project our SpatialGridDataFrame – it gets\r\ninterpolated into points, losing our nice pretty smooth heatmap.<\/p>\r\n\r\n
There are two real options for us:<\/p>\r\n\r\n