3

我试图通过覆盖更新(更高细节)卫星图像(我从 { leaflet} 包中获得)来改善 Rayshader 的外观,但覆盖与 3D 渲染不匹配。

理想情况下,我正在寻找可以获取全球卫星图像的开源解决方案。如果您找到我感兴趣的地区(夏威夷)的更详细数据,则可获得奖励积分。

geoviz使用 { } 和 { } 的一种方法使用rayshaderslippy_overlay()函数从Mapbox卫星mapbox-streets-v8mapbox-terrain-v2mapbox-traffic-v1terrain-rgbmapbox-incidents- v1 ) 或雄蕊。虽然我发现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工作流程以获取leaflet的覆盖,但结果非常奇怪。Will 的方法有点不同,因为它从 USGS 获取高程数据,它没有必须的海湾高程 - 所以我使用了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基本地图的叠加层:

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

add_overlay(., overlay_img, alpha = 0.8) 中的错误:参数 3 匹配多个形式参数

4

1 回答 1

0

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

于 2019-11-18T21:35:48.457 回答