R: Cara membuat peta panas dengan paket leaflet

10

Saya membaca posting tentang peta interaktif dengan R menggunakan leafletpaket.

Pada artikel ini, penulis membuat peta panas seperti ini:

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)

Saya tidak terbiasa dengan bkde2Dfungsinya, jadi saya bertanya-tanya apakah kode ini dapat digeneralisasi ke sembarang bentuk?

Bagaimana jika setiap node memiliki bobot tertentu, yang ingin kami wakili pada peta panas?

Apakah ada cara lain untuk membuat peta panas dengan leafletpeta di R?

Felipe
sumber
bke2d memungkinkan Anda melakukan binning 2d (estimasi kepadatan kernel) untuk satu set poin (sehingga pasangan lng / lat bekerja dengan baik). paket ks mendukung perataan kernel untuk data dari 1 hingga 6 dimensi. paket akima dapat melakukan interpolasi (berguna saat Anda membutuhkan kisi biasa). mungkin layak membaca pada tampilan tugas spasial untuk ini sebelum mencoba untuk menghasilkan sesuatu yang mungkin tidak mewakili data dengan baik.
hrbrmstr
ok, terima kasih untuk tautannya, saya pasti akan melihatnya. Sebenarnya fungsi bke2d tidak berfungsi dengan baik dengan data saya karena berfungsi dalam contoh, dan saya tidak dapat menemukan alasannya.
Felipe

Jawaban:

10

Inilah pendekatan saya untuk membuat peta panas yang lebih umum di Leaflet menggunakan R. Pendekatan ini menggunakan contourLines, seperti posting blog yang disebutkan sebelumnya, tetapi saya gunakan lapplyuntuk mengulangi semua hasil dan mengubahnya menjadi poligon umum. Pada contoh sebelumnya terserah kepada pengguna untuk secara individual memplot setiap poligon, jadi saya akan menyebutnya "lebih umum" (setidaknya ini generalisasi yang saya inginkan ketika saya membaca posting 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])

Inilah yang akan Anda miliki pada saat ini: masukkan deskripsi gambar di sini

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

Dan inilah tampilan peta panas dengan titik-titik:

masukkan deskripsi gambar di sini

Berikut adalah area yang menyarankan kepada saya bahwa saya perlu mencari beberapa parameter atau mungkin menggunakan kernel yang berbeda:

masukkan deskripsi gambar di sini

## 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
sumber
Menjelajahi jalinan yang mencoba mencari tahu ini dan sejauh ini adalah contoh terbaik yang saya temukan. Dicolokkan ke kode saya dan "hanya berfungsi." Luar biasa. Terima kasih!
Jeff Allen
Terima kasih! Saya sebenarnya telah membuat repo dengan beberapa contoh peta lain yang mungkin berguna untuk orang lain github.com/geneorama/wnv_map_demo
geneorama
Terima kasih untuk tutorial mini ini. Bagaimana Anda memilih bandwidthin bkde2d()?
the_darkside
2
@the_darkside pertanyaan bagus. Pada kenyataannya saya mengutak-atiknya sampai saya mendapatkan sesuatu yang saya sukai, saya awalnya mengembangkan peta ini secara khusus untuk memeriksa asumsi bandwidth. Dalam hal ini saya benar-benar menggunakan MASS::bandwidth.nrd(dat$latitude)dan MASS::bandwidth.nrd(dat$longitude)sebagai titik awal. Lihat ?MASS::kde2ddokumentasi yang terhubung dengan bandwith.nrd. Lihat juga ?KernSmooth::dpikapakah Anda tertarik dengan pendekatan lain.
geneorama
apakah gridsize = c(100,100)itu berarti ada total 10.000 sel?
the_darkside
4

Membangun dari jawaban genorama di atas, Anda juga dapat mengubah output bkde2D menjadi raster daripada garis kontur, menggunakan nilai fhat sebagai nilai sel 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")

Ini adalah output Anda. Perhatikan bahwa nilai kerapatan rendah masih muncul sebagai berwarna di raster.

Output Raster

Kami dapat menghapus sel dengan kerapatan rendah ini dengan yang berikut:

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

Sekarang setiap sel raster dengan nilai kurang dari 1 transparan.

Peta akhir

Jika Anda ingin raster yang di biner, gunakan fungsi colorBin daripada fungsi 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")

Kepadatan Kernel Raster Binned

Untuk membuatnya lebih lancar, cukup tambahkan ukuran grid pada fungsi bkde2D. Ini meningkatkan resolusi raster yang dihasilkan. (Saya mengubahnya menjadi

gridsize = c(1000,1000)

Keluaran:

Raster yang dihaluskan

Jacob Sanua
sumber
Bagaimana Anda bisa mengubah deskripsi legenda "Kernel Density of Points" menjadi sesuatu yang lebih intuitif, seperti "Pencurian per km persegi"? Saya kira ada persamaan yang menghubungkan bandwidth, ukuran grid dan proyeksi, atau mungkin bahkan $ khat yang menggambarkan unit.
kelima
3

Cara mudah membuat peta panas Leaflet di R adalah menggunakan plugin Leaflet.heat . Panduan yang sangat baik tentang cara menggunakannya dapat ditemukan di sini . Semoga bermanfaat.

Sam
sumber