R: Como criar um mapa de calor com o pacote do folheto

10

Eu li um post sobre mapas interativos com R usando o leafletpacote.

Neste artigo, o autor cria um mapa de calor como este:

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)

Eu não estou familiarizado com a bkde2Dfunção, então eu estou querendo saber se esse código pode ser generalizado para qualquer shapefiles?

E se cada nó tiver um peso específico, que gostaríamos de representar no mapa de calor?

Existem outras maneiras de criar um mapa de calor com o leafletmapa em R?

Felipe
fonte
O bke2d permite fazer bin 2d (estimativa da densidade do kernel) para um conjunto de pontos (para que os pares lng / lat funcionem bem). o pacote ks suporta suavização do kernel para dados de 1 a 6 dimensões. o pacote akima pode fazer interpolação (útil quando você precisa de uma grade regular). pode valer a pena ler sobre a visualização da tarefa espacial antes de tentar produzir algo que pode não representar os dados corretamente.
Hrbrmstr 03/11/2015
ok, obrigado pelo link, eu definitivamente vou procurar isso. Na verdade, a função bke2d não está funcionando tão bem com meus dados como funciona no exemplo, e não consigo entender o porquê.
Felipe

Respostas:

10

Aqui está minha abordagem para criar um mapa de calor mais generalizado no Leaflet usando R. Essa abordagem usa contourLines, como a postagem de blog mencionada anteriormente, mas eu uso lapplypara iterar todos os resultados e convertê-los em polígonos gerais. No exemplo anterior, cabe ao usuário plotar individualmente cada polígono, então eu chamaria isso de "mais generalizado" (pelo menos essa é a generalização que eu queria quando li a postagem do 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])

Aqui está o que você terá neste momento: insira a descrição da imagem aqui

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

E é assim que o mapa de calor com pontos seria:

insira a descrição da imagem aqui

Aqui está uma área que me sugere que eu preciso ajustar alguns parâmetros ou talvez usar um kernel diferente:

insira a descrição da imagem aqui

## 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
fonte
Vasculhou as interwebs tentando descobrir isso e este foi de longe o melhor exemplo que encontrei. Liguei-o ao meu código e "apenas funcionou". Impressionante. Obrigado!
Jeff Allen
Obrigado! Eu realmente criado um repo com vários outros exemplos de mapas que pode ser útil para outros github.com/geneorama/wnv_map_demo
geneorama
Obrigado por este mini-tutorial. Como você selecionou o bandwidthin bkde2d()?
the_darkside
2
@the_darkside ótima pergunta. Na realidade, eu brinco com isso até conseguir algo que gosto, originalmente desenvolvi este mapa especificamente para examinar as suposições de largura de banda. Neste caso, eu realmente usei MASS::bandwidth.nrd(dat$latitude)e MASS::bandwidth.nrd(dat$longitude)como ponto de partida. Consulte a ?MASS::kde2ddocumentação para a qual está vinculado bandwith.nrd. Veja também ?KernSmooth::dpikse você está interessado em outra abordagem.
geneorama
se gridsize = c(100,100)isso significa que há um total de 10.000 células?
21718
4

Com base na resposta do genorama acima, você também pode converter a saída do bkde2D em uma varredura em vez de linhas de contorno, usando os valores fhat como os valores das células 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")

Esta é a sua saída. Observe que os valores de baixa densidade ainda aparecem como coloridos na varredura.

Saída Raster

Podemos remover essas células de baixa densidade com o seguinte:

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

Agora qualquer célula de varredura com um valor menor que 1 é transparente.

Mapa final

Se você deseja uma varredura em bin, use a função colorBin em vez da função 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")

Densidade do kernel de varredura binada

Para torná-lo mais suave, basta aumentar o tamanho da grade na função bkde2D. Isso aumenta a resolução da varredura gerada. (Mudei para

gridsize = c(1000,1000)

Resultado:

Varredura suavizada

Jacob Sanua
fonte
Como você pode converter a descrição da legenda "Densidade de pontos do kernel" em algo mais intuitivo, como "Roubos por km quadrado"? Acho que existe uma equação que liga a largura de banda, o tamanho da grade e a projeção, ou talvez até kdf $ fhat que descreva as unidades.
fifthace 14/02
3

Uma maneira fácil de criar mapas de calor do Leaflet no R é usar o plug- in Leaflet.heat . Um excelente guia sobre como usá-lo pode ser encontrado aqui . Espero que você ache útil.

Sam
fonte