tylermorganwall / rayshader

R Package for 2D and 3D mapping and data visualization
https://www.rayshader.com/
2.05k stars 211 forks source link

Error in rayrender::obj_model(cache_filename, x = -bbox_center[1], y = -bbox_center[2] #270

Closed tonynick closed 1 year ago

tonynick commented 1 year ago

Describe the bug I am trying to replicate code to visualize population density of American states using the Kontur population dataset and rayshader. Everything works with the code until I try to render the image. I get the following error:

"Error in rayrender::obj_model(cache_filename, x=-bbox_center[1],y=-bbox_center[2]"

I tried a couple of different ways of running the code, and I tried it on a PC and on a Mac, and I got the same error.

Session Info R version 4.2.1 (2022-06-23 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale: [1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8 [4] LC_NUMERIC=C LC_TIME=English_United States.utf8

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

other attached packages: [1] colorspace_2.0-3 MetBrewer_0.2.0 rayrender_0.29.4 rayshader_0.24.10 stars_0.6-0 abind_1.4-5
[7] forcats_0.5.2 stringr_1.4.1 purrr_0.3.4 readr_2.1.2 tidyr_1.2.0 tibble_3.1.8
[13] ggplot2_3.4.0 tidyverse_1.3.2 tigris_2.0.1 sf_1.0-9 magick_2.7.3 rgl_1.0.1
[19] dplyr_1.0.9

loaded via a namespace (and not attached): [1] httr_1.4.4 foreach_1.5.2 jsonlite_1.8.0 modelr_0.1.9 assertthat_0.2.1
[6] googlesheets4_1.0.1 cellranger_1.1.0 progress_1.2.2 Rttf2pt1_1.3.10 pillar_1.8.1
[11] backports_1.4.1 glue_1.6.2 extrafontdb_1.0 uuid_1.1-0 digest_0.6.29
[16] rvest_1.0.3 htmltools_0.5.4 pkgconfig_2.0.3 broom_1.0.3 haven_2.5.1
[21] scales_1.2.1 tzdb_0.3.0 proxy_0.4-27 googledrive_2.0.0 farver_2.1.1
[26] generics_0.1.3 ellipsis_0.3.2 pacman_0.5.1 withr_2.5.0 cli_3.6.0
[31] magrittr_2.0.3 crayon_1.5.2 readxl_1.4.1 fs_1.5.2 fansi_1.0.3
[36] doParallel_1.0.17 lwgeom_0.2-8 xml2_1.3.3 class_7.3-20 prettyunits_1.1.1
[41] tools_4.2.1 hms_1.1.2 gargle_1.2.1 lifecycle_1.0.3 munsell_0.5.0
[46] reprex_2.0.2 compiler_4.2.1 e1071_1.7-11 rlang_1.0.6 classInt_0.4-7
[51] units_0.8-0 grid_4.2.1 iterators_1.0.14 rstudioapi_0.14 rappdirs_0.3.3
[56] htmlwidgets_1.6.1 base64enc_0.1-3 codetools_0.2-18 gtable_0.3.1 curl_4.3.2
[61] DBI_1.1.3 R6_2.5.1 lubridate_1.8.0 knitr_1.42 fastmap_1.1.0
[66] extrafont_0.18 utf8_1.2.2 KernSmooth_2.23-20 stringi_1.7.8 parallel_4.2.1
[71] Rcpp_1.0.9 png_0.1-8 vctrs_0.5.0 dbplyr_2.2.1 tidyselect_1.2.0
[76] xfun_0.36

Reproducible Example Please provide a reproducible example that triggers the bug below. It will not allow me to upload the GPKG file. It is available for download here.

library(sf)
library(tidyverse)
library(elevatr)
library(rayshader)
library(glue)
library(colorspace)
library(tigris)
library(stars)
library(NatParksPalettes)

# Set map name that will be used in file names, and 
# to get get boundaries from master NPS list

map <- "iowa"

# Kontur data source: https://data.humdata.org/organization/kontur
d_layers <- st_layers("data/kontur/kontur_population_US_20220630.gpkg")
d_crs <- d_layers[["crs"]][[1]][[2]]

s <- states() |> 
  st_transform(crs = d_crs)

st <- s |> 
  filter(NAME == str_to_title(str_replace_all("iowa", "_", " ")))

wkt_st <- st_as_text(st[[1,"geometry"]])

data <- st_read("data/kontur/kontur_population_US_20220630.gpkg",
                wkt_filter = wkt_st)

data |> 
  ggplot() +
  geom_sf()

st_d <- st_join(data, st, left = FALSE)

st_d |> 
  ggplot() +
  geom_sf()

#st_d <- st_intersection(data, st)

bb <- st_bbox(st_d)
yind <- st_distance(st_point(c(bb[["xmin"]], bb[["ymin"]])), 
                    st_point(c(bb[["xmin"]], bb[["ymax"]])))
xind <- st_distance(st_point(c(bb[["xmin"]], bb[["ymin"]])), 
                    st_point(c(bb[["xmax"]], bb[["ymin"]])))

if (yind > xind) {
  y_rat <- 1
  x_rat <- xind / yind
} else {
  x_rat <- 1
  y_rat <- yind / xind
}

size <- 9000
rast <- st_rasterize(st_d |> 
                       select(population, geom),
                     nx = floor(size * x_rat), ny = floor(size * y_rat))

mat <- matrix(rast$population, nrow = floor(size * x_rat), ncol = floor(size * y_rat))

# set up color palette

pal <- "olympic"

c1 <- natparks.pals("Olympic", n = 10)
colors <- c1[c(5:1, 10:6)]
swatchplot(colors)

texture <- grDevices::colorRampPalette(colors, bias = 4)(256)

swatchplot(texture)

###################
# Build 3D Object #
###################

# Keep this line so as you're iterating you don't forget to close the
# previous window

try(rgl::rgl.close())

# Create the initial 3D object
mat |> 
  height_shade(texture = texture) |> 
  plot_3d(heightmap = mat, 
          # This is my preference, I don't love the `solid` in most cases
          solid = FALSE,
          soliddepth = 100,
          # You might need to hone this in depending on the data resolution;
          # lower values exaggerate the height
          z = 12,
          # Set the location of the shadow, i.e. where the floor is.
          # This is on the same scale as your data, so call `zelev` to see the
          # min/max, and set it however far below min as you like.
          shadowdepth = 0,
          # Set the window size relatively small with the dimensions of our data.
          # Don't make this too big because it will just take longer to build,
          # and we're going to resize with `render_highquality()` below.
          windowsize = c(800,800), 
          # This is the azimuth, like the angle of the sun.
          # 90 degrees is directly above, 0 degrees is a profile view.
          phi = 90, 
          zoom = 1, 
          # `theta` is the rotations of the map. Keeping it at 0 will preserve
          # the standard (i.e. north is up) orientation of a plot
          theta = 0, 
          background = "white") 

# Use this to adjust the view after building the window object
render_camera(phi = 35, zoom = .75, theta = 20)

###############################
# Create High Quality Graphic #
###############################

# You should only move on if you have the object set up
# as you want it, including colors, resolution, viewing position, etc.

# Ensure dir exists for these graphics
if (!dir.exists(glue("images/{map}"))) {
  dir.create(glue("images/{map}"))
}

# Set up outfile where graphic will be saved.
# Note that I am not tracking the `images` directory, and this
# is because these files are big enough to make tracking them on
# GitHub difficult. 
outfile <- str_to_lower(glue("images/{map}/{map}_{pal}.png"))

# Now that everything is assigned, save these objects so we
# can use then in our markup script
saveRDS(list(
  map = map,
  pal = pal,
  colors = colors,
  outfile = outfile,
  coords = coords
), glue("R/portraits/{map}/header.rds"))

{
  # Test write a PNG to ensure the file path is good.
  # You don't want `render_highquality()` to fail after it's 
  # taken hours to render.
  if (!file.exists(outfile)) {
    png::writePNG(matrix(1), outfile)
  }
  # I like to track when I start the render
  start_time <- Sys.time()
  cat(glue("Start Time: {start_time}"), "\n")
  render_highquality(
    # We test-wrote to this file above, so we know it's good
    outfile, 
    # See rayrender::render_scene for more info, but best
    # sample method ('sobol') works best with values over 256
    samples = 450, 
    preview = FALSE,
    light = TRUE,
    lightdirection = rev(c(140, 140, 150, 150)),
    lightcolor = c(colors[1], "white", colors[9], "white"),
    lightintensity = c(750, 70, 1000, 70),
    lightaltitude = c(10, 80, 10, 80),
    # All it takes is accidentally interacting with a render that takes
    # hours in total to decide you NEVER want it interactive
    interactive = FALSE,
    # HDR lighting used to light the scene
    # environment_light = "assets/env/phalzer_forest_01_4k.hdr",
    # # environment_light = "assets/env/small_rural_road_4k.hdr",
    # # Adjust this value to brighten or darken lighting
    # intensity_env = 1.5,
    # # Rotate the light -- positive values move it counter-clockwise
    # rotate_env = 130,
    # This effectively sets the resolution of the final graphic,
    # because you increase the number of pixels here.
    # width = round(6000 * wr), height = round(6000 * hr),
    width = 9000, height = 9000
  )
  end_time <- Sys.time()
  cat(glue("Total time: {end_time - start_time}"), "\n")
}
tonynick commented 1 year ago

Oh and sorry, here is the source code I'm trying to replicate.

tonynick commented 1 year ago

Reinstalling your package using remotes::install_github seemed to fix it.