leaflet - 校正 Rayshader 的卫星图像覆盖

Question

我试图通过叠加更新的(更详细的)卫星图像(我从 {leaflet} 包中获取)来改善 Rayshader 的外观，但叠加层不匹配与 3D 渲染。

理想情况下，我正在寻找可以获取全局卫星图像的开源解决方案。如果您为我感兴趣的区域(夏威夷)找到更详细的数据，可加分。

使用 {geoviz} 和 {rayshader} 的一种方法是使用 slippy_overlay() 函数从 < strong>Mapbox(卫星，mapbox-streets-v8，mapbox-terrain-v2，mapbox-traffic- v1、terrain-rgb、mapbox-incidents-v1) 或 Stamen。虽然我发现 mapbox-terrain-v2 是最好的，但它仍然缺少我想要的细节。因为它需要为 mapbox 设置 API，所以我只使用下面的 stamen/watercolor:

library(geoviz)
library(rayshader)

### Maui
lat = 20.785700
lon = -156.259204
square_km = 22
max_tiles = 10

dem <- mapzen_dem(lat, lon, square_km, max_tiles)

elev_matrix  = matrix(
        raster::extract(dem, raster::extent(dem), buffer=1000),
        nrow = ncol(dem),
        ncol = nrow(dem)
)

ambmat <- ambient_shade(elev_matrix, zscale = 30)
raymat <- ray_shade(elev_matrix, zscale = 30, lambert = TRUE)
watermap <- detect_water(elev_matrix)

overlay_img <-
        slippy_overlay(dem,
                       image_source = "stamen",
                       image_type = "watercolor",
                       png_opacity = 0.3,
                       max_tiles = max_tiles)

elev_matrix %>%
        sphere_shade(sunangle = 270, texture = "imhof4") %>%
        add_water(detect_water(elev_matrix), color="imhof4") %>%
        add_shadow(ray_shade(elev_matrix,zscale=3,maxsearch = 300),0.5) %>%
        add_shadow(ambmat,0.5) %>%
        add_overlay(overlay_img) %>%
        plot_3d(elev_matrix,
                solid = T,
                water = T,
                waterdepth = 0,
                wateralpha = 0.5,
                watercolor = "lightblue",
                waterlinecolor = "white",
                waterlinealpha = 0.5,
                zscale= raster_zscale(dem) / 3,
                fov=0,theta=135,zoom=0.75,phi=45, windowsize = c(1000,800))

我正在尝试适应 Will Bishop的 workflow for getting overlays with the leaflet package但结果很奇怪。 Will 的方法有点不同，因为它从 USGS 获取海拔数据，它没有 baythmetric elevation，这对我来说是必须 - 所以我使用了 geoviz

library(leaflet)

# define bounding box with longitude/latitude coordinates
bbox <- list(
        p1 = list(long = -156.8037, lat = 20.29737),
        p2 = list(long = -155.7351, lat = 21.29577)
)

leaflet() %>%
        addTiles() %>% 
        addRectangles(
                lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
                lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
                fillColor = "transparent"
        ) %>%
        fitBounds(
                lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
                lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
        )

geoviz 中的山体阴影面积是多少？

dim(dem)
780 780   1

好的，覆盖图像需要是 780 x 780，所以我修改了辅助函数以下载带有 World_Imagery basemap 的覆盖:

define_image_size <- function(bbox, major_dim = 780) {
        # calculate aspect ration (width/height) from lat/long bounding box
        aspect_ratio <- abs((bbox$p1$long - bbox$p2$long) / (bbox$p1$lat - bbox$p2$lat))
        # define dimensions
        img_width <- ifelse(aspect_ratio > 1, major_dim, major_dim*aspect_ratio) %>% round()
        img_height <- ifelse(aspect_ratio < 1, major_dim, major_dim/aspect_ratio) %>% round()
        size_str <- paste(img_width, img_height, sep = ",")
        list(height = img_height, width = img_width, size = size_str)
}

get_arcgis_map_image <- function(bbox, map_type = "World_Imagery", file = NULL, 
                                 width = 780, height = 780, sr_bbox = 4326) {
        require(httr)
        require(glue) 
        require(jsonlite)

        url <- parse_url("https://utility.arcgisonline.com/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute")

        # define JSON query parameter
        web_map_param <- list(
                baseMap = list(
                        baseMapLayers = list(
                                list(url = jsonlite::unbox(glue("https://services.arcgisonline.com/ArcGIS/rest/services/{map_type}/MapServer",
                                                                map_type = map_type)))
                        )
                ),
                exportOptions = list(
                        outputSize = c(width, height)
                ),
                mapOptions = list(
                        extent = list(
                                spatialReference = list(wkid = jsonlite::unbox(sr_bbox)),
                                xmax = jsonlite::unbox(max(bbox$p1$long, bbox$p2$long)),
                                xmin = jsonlite::unbox(min(bbox$p1$long, bbox$p2$long)),
                                ymax = jsonlite::unbox(max(bbox$p1$lat, bbox$p2$lat)),
                                ymin = jsonlite::unbox(min(bbox$p1$lat, bbox$p2$lat))
                        )
                )
        )

        res <- GET(
                url, 
                query = list(
                        f = "json",
                        Format = "PNG32",
                        Layout_Template = "MAP_ONLY",
                        Web_Map_as_JSON = jsonlite::toJSON(web_map_param))
        )

        if (status_code(res) == 200) {
                body <- content(res, type = "application/json")
                message(jsonlite::toJSON(body, auto_unbox = TRUE, pretty = TRUE))
                if (is.null(file)) 
                        file <- tempfile("overlay_img", fileext = ".png")

                img_res <- GET(body$results[[1]]$value$url)
                img_bin <- content(img_res, "raw")
                writeBin(img_bin, file)
                message(paste("image saved to file:", file))
        } else {
                message(res)
        }
        invisible(file)
}

现在下载文件，然后加载它

image_size <- define_image_size(bbox, major_dim = 780)

# fetch overlay image
overlay_file <- "maui_overlay.png"
get_arcgis_map_image(bbox, map_type = "World_Imagery", file = overlay_file,
                     # width = image_size$width, height = image_size$height, 
                     sr_bbox = 4326)

overlay_img <- png::readPNG("maui_overlay.png")

好吧，让我们开始吧

elev_matrix %>%
        sphere_shade(sunangle = 270, texture = "imhof4") %>%
        add_water(detect_water(elev_matrix), color="imhof4") %>%
        add_shadow(ray_shade(elev_matrix,zscale=3,maxsearch = 300),0.5) %>%
        add_shadow(ambmat,0.5) %>%
        add_overlay(overlay_img, alphacolor = 1) %>%
        plot_3d(elev_matrix,
                solid = T,
                water = T,
                waterdepth = 0,
                wateralpha = 0.5,
                watercolor = "lightblue",
                waterlinecolor = "white",
                waterlinealpha = 0.5,
                zscale= raster_zscale(dem) / 3,
                fov=0,theta=135,zoom=0.75,phi=45, windowsize = c(1000,800))

如您所见，叠加图像已旋转到山体阴影。

现在我也意识到，当您尝试显示深海矩阵数据时，使用边界框方法获取卫星并不理想。以某种方式以编程方式对该叠加层进行子集化是理想的，但一旦我弄清楚如何旋转叠加层，我可能最终会使用 inkscape。

我尝试使用 {magick} 的 image_rotate() 函数无济于事:

library(magick)
maui <- magick::image_read("maui_overlay.png")
image_rotate(maui, 30) # -> maui_30
# image_write(maui_30, path = "maui_overlay_30.png", format = "png")

但是 magick 改变了维度:

# A tibble: 1 x 7
  format width height colorspace matte filesize density
  <chr>  <int>  <int> <chr>      <lgl>    <int> <chr>  
1 PNG     1068   1068 sRGB       TRUE         0 38x38  

rayshader 会报错:

overlay_img <- png::readPNG("maui_overlay_30.png")
elev_matrix %>%
        sphere_shade(sunangle = 270, texture = "imhof4") %>%
        add_water(detect_water(elev_matrix), color="imhof4") %>%
        add_shadow(ray_shade(elev_matrix,zscale=3,maxsearch = 300),0.5) %>%
        add_shadow(ambmat,0.5) %>%
        add_overlay(overlay_img, alphacolor = 1) %>%
        plot_3d(elev_matrix,
                solid = T,
                water = T,
                waterdepth = 0,
                wateralpha = 0.5,
                watercolor = "lightblue",
                waterlinecolor = "white",
                waterlinealpha = 0.5,
                zscale= raster_zscale(dem) / 3,
                fov=0,theta=135,zoom=0.75,phi=45, windowsize = c(1000,800))

Error in add_overlay(., overlay_img, alpha = 0.8) : argument 3 matches multiple formal arguments

Answer 1

答案再简单不过了……它需要转置 overlay_img = aperm(overlay_img, c(2,1,3))。