R: Wie erstelle ich eine Heatmap mit dem Faltblattpaket?

10

Ich habe mit dem Paket einen Beitrag über interaktive Karten mit R gelesen leaflet.

In diesem Artikel hat der Autor eine Heatmap wie folgt erstellt:

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)

Ich bin mit der bkde2DFunktion nicht vertraut und frage mich, ob dieser Code auf Shapefiles verallgemeinert werden kann.

Was ist, wenn jeder Knoten ein bestimmtes Gewicht hat, das wir auf der Heatmap darstellen möchten?

Gibt es andere Möglichkeiten, eine Heatmap mit leafletMap in R zu erstellen ?

Felipe
quelle
Mit bke2d können Sie 2D-Binning (Kernel-Dichteschätzung) für eine Reihe von Punkten durchführen (lng / lat-Paare funktionieren also gut). Das ks-Paket unterstützt die Kernel-Glättung für Daten von 1 bis 6 Dimensionen. Das Akima-Paket kann interpolieren (nützlich, wenn Sie ein reguläres Raster benötigen). Es kann sich lohnen, die räumliche Aufgabenansicht zu lesen, bevor Sie versuchen, etwas zu erstellen, das die Daten möglicherweise nicht richtig darstellt.
hrbrmstr
ok, danke für den Link, ich werde das auf jeden Fall sehen. Tatsächlich funktioniert die Funktion bke2d mit meinen Daten nicht so gut wie im Beispiel, und ich kann mir nicht vorstellen, warum.
Felipe

Antworten:

10

Hier ist mein Ansatz zum Erstellen einer allgemeineren Heatmap in Leaflet mit R. Dieser Ansatz verwendet contourLineswie der zuvor erwähnte Blog-Beitrag, aber ich lapplyiteriere über alle Ergebnisse und konvertiere sie in allgemeine Polygone. Im vorherigen Beispiel ist es Sache des Benutzers, jedes Polygon einzeln zu zeichnen, daher würde ich dies als "allgemeiner" bezeichnen (zumindest ist dies die Generalisierung, die ich beim Lesen des Blogposts wollte!).

## 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])

Folgendes haben Sie an dieser Stelle: Geben Sie hier die Bildbeschreibung ein

## 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")

Und so würde die Heatmap mit Punkten aussehen:

Geben Sie hier die Bildbeschreibung ein

Hier ist ein Bereich, der mir vorschlägt, dass ich einige Parameter optimieren oder vielleicht einen anderen Kernel verwenden muss:

Geben Sie hier die Bildbeschreibung ein

## 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])
Genorama
quelle
Durchsuchte die Interwebs, um dies herauszufinden, und dies war bei weitem das beste Beispiel, das ich gefunden habe. Steckte es in meinen Code und es "funktionierte einfach". Genial. Vielen Dank!
Jeff Allen
Vielen Dank! Ich habe tatsächlich ein Repo mit mehreren anderen Karten Beispiele erstellt , die für andere nützlich sein könnten github.com/geneorama/wnv_map_demo
geneorama
Danke für dieses Mini-Tutorial. Wie haben Sie das bandwidthin ausgewählt bkde2d()?
the_darkside
2
@the_darkside tolle Frage. In Wirklichkeit spiele ich damit herum, bis ich etwas bekomme, das mir gefällt. Ich habe diese Karte ursprünglich speziell entwickelt, um die Bandbreitenannahmen zu untersuchen. In diesem Fall habe ich tatsächlich MASS::bandwidth.nrd(dat$latitude)und MASS::bandwidth.nrd(dat$longitude)als Ausgangspunkt verwendet. Siehe ?MASS::kde2dDokumentation, auf die verwiesen wird bandwith.nrd. Sehen ?KernSmooth::dpikSie auch nach, ob Sie an einem anderen Ansatz interessiert sind.
Genorama
Wenn das gridsize = c(100,100)bedeutet, dass es insgesamt 10.000 Zellen gibt?
the_darkside
4

Aufbauend auf der obigen Antwort von genorama können Sie auch die Ausgabe von bkde2D in ein Raster anstatt in Konturlinien konvertieren, wobei Sie die fhat-Werte als Rasterzellenwerte verwenden

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")

Dies ist Ihre Ausgabe. Beachten Sie, dass die Werte mit niedriger Dichte im Raster immer noch farbig angezeigt werden.

Rasterausgabe

Wir können diese Zellen mit niedriger Dichte folgendermaßen entfernen:

#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")

Jetzt ist jede Rasterzelle mit einem Wert von weniger als 1 transparent.

Letzte Karte

Wenn Sie ein gruppiertes Raster möchten, verwenden Sie die Funktion colorBin anstelle der Funktion 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")

Binned Raster Kernel Density

Um es glatter zu machen, erhöhen Sie einfach die Rastergröße in der bkde2D-Funktion. Dies erhöht die Auflösung des generierten Rasters. (Ich habe es geändert in

gridsize = c(1000,1000)

Ausgabe:

Glättungsraster

Jacob Sanua
quelle
Wie können Sie die Legendenbeschreibung „Kernel Density of Points“ in etwas Intuitiveres wie „Diebstähle pro Quadratkilometer“ umwandeln? Ich denke, es gibt eine Gleichung, die Bandbreite, Rastergröße und Projektion verbindet, oder vielleicht sogar kdf $ fhat, das die Einheiten beschreibt.
Fifthace
3

Eine einfache Möglichkeit, Leaflet- Heatmaps in R zu erstellen, ist das Leaflet.heat- Plugin. Eine ausgezeichnete Anleitung zur Verwendung finden Sie hier . Ich hoffe, Sie finden es nützlich.

Sam
quelle