使用 stat_density_2d(geom = "polygon") 时没有渲染

Nothing rendering when using stat_density_2d(geom = "polygon")

出于某种原因,在指定 geom = "polygon" 时,stat_density_2d() 似乎不适合我,我完全被难住了。这是我的代码...

library(sf)
library(tidyverse)
library(RANN2)
library(hexbin)
library(mapproj)

options(stringsAsFactors = FALSE)

raleigh_police <- rgdal::readOGR("https://opendata.arcgis.com/datasets/24c0b37fa9bb4e16ba8bcaa7e806c615_0.geojson", "OGRGeoJSON")

raleigh_police_sf <- raleigh_police %>%
  st_as_sf()

raleigh_police_sf %>%
    filter(crime_description == "Burglary/Residential") %>%
    st_coordinates() %>%
    as_tibble() %>%
    ggplot() +
    stat_density_2d(aes(X, Y, fill = stat(level)), geom = "polygon")

这是我的 sessionInfo()...

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.1

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2.2  mapproj_1.2.6   maps_3.3.0      hexbin_1.27.2   RANN2_0.1       forcats_0.3.0   stringr_1.3.1   dplyr_0.7.8    
 [9] purrr_0.2.5     readr_1.1.1     tidyr_0.8.2     tibble_1.4.2    ggplot2_3.1.0   tidyverse_1.2.1 sf_0.7-1       

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0       lubridate_1.7.4  lattice_0.20-38  class_7.3-14     utf8_1.1.4       assertthat_0.2.0 rprojroot_1.3-2 
 [8] digest_0.6.18    R6_2.3.0         cellranger_1.1.0 plyr_1.8.4       backports_1.1.2  evaluate_0.12    e1071_1.7-0     
[15] httr_1.3.1       blogdown_0.9     pillar_1.3.0     rlang_0.3.0.1    lazyeval_0.2.1   readxl_1.1.0     rstudioapi_0.8  
[22] rmarkdown_1.10   labeling_0.3     rgdal_1.3-6      munsell_0.5.0    broom_0.5.0      compiler_3.5.1   modelr_0.1.2    
[29] xfun_0.4         pkgconfig_2.0.2  htmltools_0.3.6  tidyselect_0.2.5 bookdown_0.7     codetools_0.2-15 fansi_0.4.0     
[36] crayon_1.3.4     withr_2.1.2      MASS_7.3-51.1    grid_3.5.1       nlme_3.1-137     spData_0.2.9.4   jsonlite_1.5    
[43] gtable_0.2.0     DBI_1.0.0        magrittr_1.5     units_0.6-1      scales_1.0.0     cli_1.0.1        stringi_1.2.4   
[50] sp_1.3-1         xml2_1.2.0       tools_3.5.1      glue_1.3.0       hms_0.4.2        yaml_2.2.0       colorspace_1.3-2
[57] classInt_0.2-3   rvest_0.3.2      knitr_1.20       bindr_0.1.1      haven_1.1.2 

我只是得到一个空白图,上面什么也没有。完全难住了。我在这里做错了什么?

更新 (2018-11-18)

原来这里的主要问题是options(stringsAsFactors = FALSE)。如果您将其注释掉并 运行 原始代码,实际上一切正常。我发现这个 GitHub Issue 这就是我尝试这样做的原因。这个问题的答案中提供了更有效的代码解决方案,他们还确保不使用 options(stringsAsFactors = FALSE).

如前所述,这只是点级数据,他们提供的 CSV 可以很好地替代:

library(tidyverse)

rp_csv_url <- "https://opendata.arcgis.com/datasets/24c0b37fa9bb4e16ba8bcaa7e806c615_0.csv"

httr::GET(
  url = rp_csv_url,
  httr::write_disk(basename(rp_csv_url)), # won't overwrite if it exists unless explicitly told to so you get caching for free
  httr::progress() # I suspect this is a big file so it's nice to see a progress bar
)

raleigh_police <- read_csv(basename(rp_csv_url))

mutate(
  raleigh_police, 
  longitude = as.numeric(longitude), # they come in wonky, still
  latitude = as.numeric(latitude) # they come in wonky, still
) -> raleigh_police

raleigh_police %>%
  filter(crime_description == "Burglary/Residential") %>% 
  ggplot() +
  stat_density_2d(
    aes(longitude, latitude, fill = stat(level)), 
    color = "#2b2b2b", size=0.125, geom = "polygon"
  ) +
  viridis::scale_fill_viridis(direction=-1, option="magma") +
  hrbrthemes::theme_ipsum_rc()

如果您想将 level 变成更有意义的东西:

h <- c(MASS::bandwidth.nrd(rp_br$longitude), 
       MASS::bandwidth.nrd(rp_br$latitude))

dens <- MASS::kde2d(
  rp_br$longitude, rp_br$latitude, h = h, n = 100
)

breaks <- pretty(range(dens$z), 10)

zdf <- data.frame(expand.grid(x = dens$x, y = dens$y), z = as.vector(dens$z))

z <- tapply(zdf$z, zdf[c("x", "y")], identity)

cl <- grDevices::contourLines(
  x = sort(unique(dens$x)), y = sort(unique(dens$y)), z = dens$z,
  levels = breaks
)

sp::SpatialPolygons(
  lapply(1:length(cl), function(idx) {
    sp::Polygons(
      srl = list(sp::Polygon(
        matrix(c(cl[[idx]]$x, cl[[idx]]$y), nrow=length(cl[[idx]]$x), byrow=FALSE)
      )),
      ID = idx
    )
  })
) -> cont

sp::coordinates(rp_br) <- ~longitude+latitude

然后:

data_frame(
  ct = sapply(sp::over(cont, sp::geometry(rp_br), returnList = TRUE), length),
  id = 1:length(ct),
  lvl = sapply(cl, function(x) x$level)
) %>% 
  count(lvl, wt=ct) %>% 
  mutate(
    pct = n/nrow(rp_br),
    pct_lab = sprintf("%s of the points fall within this level", scales::percent(pct))
  )
## # A tibble: 10 x 4
##      lvl     n     pct pct_lab                                   
##    <dbl> <int>   <dbl> <chr>                                     
##  1   10.  7302 0.927   92.7% of the points fall within this level
##  2   20.  6243 0.792   79.2% of the points fall within this level
##  3   30.  4786 0.607   60.7% of the points fall within this level
##  4   40.  3204 0.407   40.7% of the points fall within this level
##  5   50.  1945 0.247   24.7% of the points fall within this level
##  6   60.  1277 0.162   16.2% of the points fall within this level
##  7   70.   793 0.101   10.1% of the points fall within this level
##  8   80.   474 0.0601  6.0% of the points fall within this level 
##  9   90.   279 0.0354  3.5% of the points fall within this level 
## 10  100.    44 0.00558 0.6% of the points fall within this level 

除了在阅读和将 readOGR 更改为 read_sf 之前下载文件外,它对我来说按原样工作,保存由空几何引起的几个 NA 点的警告:

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3

path <- "~/Downloads/raleigh.geojson"
download.file(
    "https://opendata.arcgis.com/datasets/24c0b37fa9bb4e16ba8bcaa7e806c615_0.geojson",
    path,
    method = "curl"
)

raleigh_police <- sf::read_sf(path, "OGRGeoJSON")

raleigh_police %>%
    filter(crime_description == "Burglary/Residential") %>%
    st_coordinates() %>%
    as_tibble() %>%
    ggplot() +
    stat_density_2d(aes(X, Y, fill = stat(level)), geom = "polygon")
#> Warning: Removed 5 rows containing non-finite values (stat_density2d).

空行:

raleigh_police %>% 
    filter(crime_description == "Burglary/Residential", 
           st_is_empty(.))
#> Simple feature collection with 5 features and 21 fields (with 5 geometries empty)
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#> # A tibble: 5 x 22
#>   OBJECTID GlobalID case_number crime_category crime_code crime_descripti…
#>      <int> <chr>    <chr>       <chr>          <chr>      <chr>           
#> 1   205318 8057315… P14076062   BURGLARY/RESI… 30B        Burglary/Reside…
#> 2   417488 70afb27… P15027702   BURGLARY/RESI… 30B        Burglary/Reside…
#> 3   424718 bdf69fa… P18029113   BURGLARY/RESI… 30B        Burglary/Reside…
#> 4   436550 711c05b… P18044139   BURGLARY/RESI… 30B        Burglary/Reside…
#> 5   442091 9d7a008… P18051764   BURGLARY/RESI… 30B        Burglary/Reside…
#> # … with 16 more variables: crime_type <chr>, reported_block_address <chr>,
#> #   city_of_incident <chr>, city <chr>, district <chr>, reported_date <dttm>,
#> #   reported_year <int>, reported_month <int>, reported_day <int>,
#> #   reported_hour <int>, reported_dayofwk <chr>, latitude <dbl>,
#> #   longitude <dbl>, agency <chr>, updated_date <dttm>, geometry <POINT [°]>

sessionInfo()
#> R version 3.5.1 (2018-07-02)
#> Platform: x86_64-apple-darwin15.6.0 (64-bit)
#> Running under: macOS  10.14.1
#> 
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] sf_0.7-1           forcats_0.3.0      stringr_1.3.1     
#>  [4] dplyr_0.7.99.9000  purrr_0.2.5        readr_1.2.0       
#>  [7] tidyr_0.8.2        tibble_1.4.99.9005 ggplot2_3.1.0     
#> [10] tidyverse_1.2.1   
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_0.2.5  haven_1.1.2       lattice_0.20-35  
#>  [4] colorspace_1.3-2  htmltools_0.3.6   yaml_2.2.0       
#>  [7] rlang_0.3.0.1     e1071_1.7-0       pillar_1.3.0.9001
#> [10] glue_1.3.0        withr_2.1.2       DBI_1.0.0        
#> [13] modelr_0.1.2      readxl_1.1.0      plyr_1.8.4       
#> [16] munsell_0.5.0     gtable_0.2.0      cellranger_1.1.0 
#> [19] rvest_0.3.2       evaluate_0.12     labeling_0.3     
#> [22] knitr_1.20        class_7.3-14      broom_0.5.0      
#> [25] Rcpp_0.12.19.3    scales_1.0.0      backports_1.1.2  
#> [28] classInt_0.2-3    jsonlite_1.5      hms_0.4.2.9001   
#> [31] digest_0.6.18     stringi_1.2.4     grid_3.5.1       
#> [34] rprojroot_1.3-2   cli_1.0.1         tools_3.5.1      
#> [37] magrittr_1.5      lazyeval_0.2.1    crayon_1.3.4     
#> [40] pkgconfig_2.0.2   MASS_7.3-51       xml2_1.2.0       
#> [43] spData_0.2.9.4    lubridate_1.7.4   assertthat_0.2.0 
#> [46] rmarkdown_1.10    httr_1.3.1        R6_2.3.0         
#> [49] units_0.6-1       nlme_3.1-137      compiler_3.5.1