R: Comment construire une carte thermique avec le package de la brochure

10

J'ai lu un article sur les cartes interactives avec R en utilisant le leafletpackage.

Dans cet article, l'auteur a créé une carte thermique comme celle-ci:

X=cbind(lng,lat)
kde2d <- bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])))

x=kde2d$x1
y=kde2d$x2
z=kde2d$fhat
CL=contourLines(x , y , z)

m = leaflet() %>% addTiles() 
m %>% addPolygons(CL[[5]]$x,CL[[5]]$y,fillColor = "red", stroke = FALSE)

Je ne connais pas la bkde2Dfonction, je me demande donc si ce code pourrait être généralisé à des fichiers de formes?

Et si chaque nœud avait un poids spécifique, que nous aimerions représenter sur la carte thermique?

Y a-t-il d'autres façons de créer une carte thermique avec une leafletcarte en R?

Felipe
la source
bke2d vous permet de faire un binning 2d (estimation de la densité du noyau) pour un ensemble de points (donc les paires lng / lat fonctionnent bien). le package ks prend en charge le lissage du noyau pour les données de 1 à 6 dimensions. le package akima peut effectuer une interpolation (utile lorsque vous avez besoin d'une grille régulière). il peut être utile de lire la vue des tâches spatiales pour cela avant d'essayer de produire quelque chose qui pourrait ne pas représenter correctement les données.
hrbrmstr
ok, merci pour le lien, je vais certainement regarder ça. En fait, la fonction bke2d ne fonctionne pas aussi bien avec mes données que dans l'exemple, et je ne peux pas comprendre pourquoi.
Felipe

Réponses:

10

Voici mon approche pour faire une carte de chaleur plus généralisée dans Leaflet en utilisant R. Cette approche utilise contourLines, comme le blog mentionné précédemment, mais j'utilise lapplypour parcourir tous les résultats et les convertir en polygones généraux. Dans l'exemple précédent, c'est à l'utilisateur de tracer individuellement chaque polygone, donc j'appellerais cela "plus généralisé" (du moins c'est la généralisation que je voulais quand j'ai lu le billet de blog!).

## INITIALIZE
library("leaflet")
library("data.table")
library("sp")
library("rgdal")
# library("maptools")
library("KernSmooth")

inurl <- "https://data.cityofchicago.org/api/views/22s8-eq8h/rows.csv?accessType=DOWNLOAD"
infile <- "mvthefts.csv"

## LOAD DATA
## Also, clean up variable names, and convert dates
if(!file.exists(infile)){
    download.file(url = inurl, destfile = infile)
}
dat <- data.table::fread(infile)
setnames(dat, tolower(colnames(dat)))
setnames(dat, gsub(" ", "_", colnames(dat)))
dat <- dat[!is.na(longitude)]
dat[ , date := as.IDate(date, "%m/%d/%Y")]

## MAKE CONTOUR LINES
## Note, bandwidth choice is based on MASS::bandwidth.nrd()
kde <- bkde2D(dat[ , list(longitude, latitude)],
              bandwidth=c(.0045, .0068), gridsize = c(100,100))
CL <- contourLines(kde$x1 , kde$x2 , kde$fhat)

## EXTRACT CONTOUR LINE LEVELS
LEVS <- as.factor(sapply(CL, `[[`, "level"))
NLEV <- length(levels(LEVS))

## CONVERT CONTOUR LINES TO POLYGONS
pgons <- lapply(1:length(CL), function(i)
    Polygons(list(Polygon(cbind(CL[[i]]$x, CL[[i]]$y))), ID=i))
spgons = SpatialPolygons(pgons)

## Leaflet map with polygons
leaflet(spgons) %>% addTiles() %>% 
    addPolygons(color = heat.colors(NLEV, NULL)[LEVS])

Voici ce que vous aurez à ce stade: entrez la description de l'image ici

## Leaflet map with points and polygons
## Note, this shows some problems with the KDE, in my opinion...
## For example there seems to be a hot spot at the intersection of Mayfield and
## Fillmore, but it's not getting picked up.  Maybe a smaller bw is a good idea?

leaflet(spgons) %>% addTiles() %>%
    addPolygons(color = heat.colors(NLEV, NULL)[LEVS]) %>%
    addCircles(lng = dat$longitude, lat = dat$latitude,
               radius = .5, opacity = .2, col = "blue")

Et voici à quoi ressemblerait la carte de chaleur avec des points:

entrez la description de l'image ici

Voici un domaine qui me suggère que je dois régler certains paramètres ou peut-être utiliser un noyau différent:

entrez la description de l'image ici

## Leaflet map with polygons, using Spatial Data Frame
## Initially I thought that the data frame structure was necessary
## This seems to give the same results, but maybe there are some 
## advantages to using the data.frame, e.g. for adding more columns
spgonsdf = SpatialPolygonsDataFrame(Sr = spgons,
                                    data = data.frame(level = LEVS),
                                    match.ID = TRUE)
leaflet() %>% addTiles() %>%
    addPolygons(data = spgonsdf,
                color = heat.colors(NLEV, NULL)[spgonsdf@data$level])
geneorama
la source
Parcouru les interwebs en essayant de comprendre cela et ce fut de loin le meilleur exemple que j'ai trouvé. Je l'ai branché à mon code et ça "a juste fonctionné". Impressionnant. Je vous remercie!
Jeff Allen
Merci! J'ai en fait créé un dépôt avec plusieurs autres exemples de cartes qui pourraient être utiles pour d'autres github.com/geneorama/wnv_map_demo
geneorama
Merci pour ce mini-tutoriel. Comment avez - vous choisi le bandwidthdans bkde2d()?
the_darkside
2
@the_darkside grande question. En réalité, je joue avec jusqu'à ce que j'obtienne quelque chose que j'aime, j'ai initialement développé cette carte spécifiquement pour examiner les hypothèses de bande passante. Dans ce cas, j'ai utilisé MASS::bandwidth.nrd(dat$latitude)et MASS::bandwidth.nrd(dat$longitude)comme point de départ. Voir la ?MASS::kde2ddocumentation qui renvoie à bandwith.nrd. Vérifiez également ?KernSmooth::dpiksi vous êtes intéressé par une autre approche.
geneorama
si gridsize = c(100,100)cela signifie-t-il qu'il y a au total 10 000 cellules?
the_darkside
4

En s'appuyant sur la réponse de genorama ci-dessus, vous pouvez également convertir la sortie de bkde2D en raster plutôt qu'en lignes de contour, en utilisant les valeurs fhat comme valeurs de cellule raster

library("leaflet")
library("data.table")
library("sp")
library("rgdal")
# library("maptools")
library("KernSmooth")
library("raster")

inurl <- "https://data.cityofchicago.org/api/views/22s8-eq8h/rows.csv?accessType=DOWNLOAD"
infile <- "mvthefts.csv"

## LOAD DATA
## Also, clean up variable names, and convert dates
if(!file.exists(infile)){
  download.file(url = inurl, destfile = infile)
}
dat <- data.table::fread(infile)
setnames(dat, tolower(colnames(dat)))
setnames(dat, gsub(" ", "_", colnames(dat)))
dat <- dat[!is.na(longitude)]
dat[ , date := as.IDate(date, "%m/%d/%Y")]

## Create kernel density output
kde <- bkde2D(dat[ , list(longitude, latitude)],
              bandwidth=c(.0045, .0068), gridsize = c(100,100))
# Create Raster from Kernel Density output
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))

#create pal function for coloring the raster
palRaster <- colorNumeric("Spectral", domain = KernelDensityRaster@data@values)

## Leaflet map with raster
leaflet() %>% addTiles() %>% 
  addRasterImage(KernelDensityRaster, 
                 colors = palRaster, 
                 opacity = .8) %>%
  addLegend(pal = palRaster, 
            values = KernelDensityRaster@data@values, 
            title = "Kernel Density of Points")

Ceci est votre résultat. Notez que les valeurs de faible densité apparaissent toujours comme colorées dans le raster.

Sortie raster

Nous pouvons supprimer ces cellules de faible densité avec les éléments suivants:

#set low density cells as NA so we can make them transparent with the colorNumeric function
 KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA

#create pal function for coloring the raster
palRaster <- colorNumeric("Spectral", domain = KernelDensityRaster@data@values, na.color = "transparent")

## Redraw the map
leaflet() %>% addTiles() %>% 
  addRasterImage(KernelDensityRaster, 
                 colors = palRaster, 
                 opacity = .8) %>%
  addLegend(pal = palRaster, 
            values = KernelDensityRaster@data@values, 
            title = "Kernel Density of Points")

Désormais, toute cellule raster d'une valeur inférieure à 1 est transparente.

Carte finale

Si vous souhaitez un raster groupé, utilisez la fonction colorBin plutôt que la fonction colorNumeric:

palRaster <- colorBin("Spectral", bins = 7, domain = KernelDensityRaster@data@values, na.color = "transparent")

## Leaflet map with raster
leaflet() %>% addTiles() %>% 
  addRasterImage(KernelDensityRaster, 
                 colors = palRaster, 
                 opacity = .8) %>%
  addLegend(pal = palRaster, 
            values = KernelDensityRaster@data@values, 
            title = "Kernel Density of Points")

Densité du noyau raster en bac

Pour le rendre plus fluide, augmentez simplement la taille de la grille dans la fonction bkde2D. Cela augmente la résolution du raster généré. (Je l'ai changé en

gridsize = c(1000,1000)

Production:

Raster lissé

Jacob Sanua
la source
Comment pouvez-vous convertir la description de la légende «Densité du noyau de points» en quelque chose de plus intuitif, comme «Vols au km2»? Je suppose qu'il existe une équation reliant la bande passante, la taille de la grille et la projection, ou peut-être même kdf $ fhat qui décrit les unités.
fifthace
3

Un moyen simple de créer des cartes de chaleur Leaflet dans R est d'utiliser le plugin Leaflet.heat . Un excellent guide sur la façon de l'utiliser peut être trouvé ici . J'espère que vous le trouverez utile.

Sam
la source