Support for Bing maps?
Robinlovelace opened this issue · 10 comments
Microsoft's Bing maps have the best data for some places it seems, including my hometown of Hereford. I had a little play but cannot figure out how to get these tiles with ceramic. Would it be a lot of work? May be up for giving it a bash at some point...
You can definitely get the tiles, as shown in the reprex below using rosm.
# install.packages("rosm")
library(rosm)
library(prettymapr)
bb = searchbbox("hereford cathedral")
bmaps.plot(bb, zoomin = 1)
# install.packages("rosm")
library(rosm)
library(prettymapr)
bb = searchbbox("hereford cathedral")
#> Using default API key for pickpoint.io. If batch geocoding, please get your own (free) API key at https://app.pickpoint.io/sign-up
bmaps.plot(bb, zoomin = 1)
#> Copyright © 2021 Microsoft and its suppliers. All rights reserved. This API cannot be accessed and the content and any results may not be used, reproduced or transmitted in any manner without express written permission from Microsoft Corporation.
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definition
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definition
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definition
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definition
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definition
#> Zoom: 20
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definition
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definition
#> Fetching 24 missing tiles
#> | | | 0% | |=== | 4% | |====== | 8% | |========= | 12% | |============ | 17% | |=============== | 21% | |================== | 25% | |==================== | 29% | |======================= | 33% | |========================== | 38% | |============================= | 42% | |================================ | 46% | |=================================== | 50% | |====================================== | 54% | |========================================= | 58% | |============================================ | 62% | |=============================================== | 67% | |================================================== | 71% | |==================================================== | 75% | |======================================================= | 79% | |========================================================== | 83% | |============================================================= | 88% | |================================================================ | 92% | |=================================================================== | 96% | |======================================================================| 100%
#> ...complete!
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definition
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
#> +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
#> +wktext +no_defs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum WGS_1984 in Proj4 definitionCreated on 2021-10-14 by the reprex package (v2.0.1)
hey, unless you really want the tiles downloaded (that's what ceramic should really limit itself to) ... I would ignore ceramic and go with gdalio:
ex <- c(-2.7171221, -2.7146559, 52.0541044, 52.0549124)
#ex <- c(t(bb))
library(gdalio)
## set target grid (can be anything, any projection), it's just 6 numbers and a string
## can also use stars/terra/raster/grd objects
g <- gdalio_set_default_grid(list(extent = ex, dimension = c(512, 256), projection = "OGC:CRS84"))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
ve_src <- '<GDAL_WMS>
<Service name="VirtualEarth">
<ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
</Service>
<MaxConnections>4</MaxConnections>
<Cache/>
</GDAL_WMS>'
## get image (my _terra _stars _raster wrappers need updating hence band_output_type for now)
im <- gdalio_terra(ve_src, bands = 1:3, band_output_type = "integer")
terra::plotRGB(im)
Created on 2021-10-15 by the reprex package (v2.0.1)
gdal caches the tiles itself in its own format, so you'll likely get a local /gdalwmscache/ and you can configure where that goes (but that's another topic) - but this is the fastest and most flexible way of getting imagery for viewing in R
That is amazing! Many thanks... Will have a play but pretty sure this fixes the issue. Follow-up question: if there were to be any wrapper code to make this easier, fewer lines of code, where should it live? So many pkgs!
I dunno ... we need a list of gdal sources, I've tried doing that but it's not really my scene I just keep rough lists.
vapour_warp_raster() is the workhorse, you need extent, dimension, projection that's all.
gdalio is an experiment in "set my grid", but it's kind of unnecessary and is a throwaway package- it just makes hosting wrappers gdalio_raster() gdalio_stars() etc easy. For me the task is very straightforward and vapour does all the work - the rest is reformatting, and I guess gdalio is a kind of proto example of what could be done. It's about enough for me so I'm unlikely to push it further, the format-specific stuff belongs in format-specific packages, but those big player authors won't play ball with generic tools so we're stuck. I'm not going to fight that, there's no leadership or direction whatsoever but maybe that's a good thing. Nature red in tooth and claw, see who exists in a few years.
The real work for me is getting the raster logic separated out and paleolimbot/grd is really helpful for that. Getting imagery is essentially solved for me, and I'm going to focus on the mesh stuff again. Best wishes ;)
Many thanks for the follow-up comments @mdsumner all this is interesting and I'm sure the collective open source community will find solutions through communications, including this one, hopefully a small step towards what is useful. For context, this issue actually comes from my dad David, you know each other from the Manifold days many years ago, he remembers you!
Many thanks for dad. For completeness, here is a full reprex that gets dad what he needed for a project on historic landscape and biodiversity mapping in Herefordshire. The Venn diagram overlap of old school naturalists and programmers is pretty small as you can imagine so very few people can get this kind of data let alone use it in the county!
Any comments on the rough and ready code below welcome, and any further comms on the bigger picture chat welcome, thanks again : )
remotes::install_github("hypertidy/gdalio")
# ex <- c(-2.7171221, -2.7146559, 52.0541044, 52.0549124)
# # get bb interactively
# m = mapedit::editMap()
# mapview::mapview(m)
# bb = sf::st_bbox(m)
# ex = as.numeric(bb)[c(1, 3, 2, 4)]
# dput(ex)
ex = c(-2.518027, -2.431069, 52.158692, 52.188771)
m1 = matrix(ex, ncol = 2)
m2 = rbind(m1, m1[1, ])
m = sf::st_polygon(list(m2))
m = sf::st_sfc(m, crs = 4326)
m_grid = sf::st_make_grid(x = m, n = c(3, 2))
#ex <- c(t(bb))
library(gdalio)
## set target grid (can be anything, any projection), it's just 6 numbers and a string
## can also use stars/terra/raster/grd objects
g <- gdalio_set_default_grid(list(extent = ex, dimension = c(512, 256), projection = "OGC:CRS84"))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
ve_src <- '<GDAL_WMS>
<Service name="VirtualEarth">
<ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
</Service>
<MaxConnections>4</MaxConnections>
<Cache/>
</GDAL_WMS>'
## get image (my _terra _stars _raster wrappers need updating hence band_output_type for now)
im <- gdalio_terra(ve_src, bands = 1:3, band_output_type = "integer")
class(im)
terra::plotRGB(im)
terra::writeRaster(im, filename = "/tmp/cathedral.tif")
get_im_bing = function(x) {
bb = sf::st_bbox(x)
ex = as.numeric(bb)[c(1, 3, 2, 4)]
g = gdalio_set_default_grid(list(extent = ex, dimension = c(512, 256), projection = "OGC:CRS84"))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
ve_src <- '<GDAL_WMS>
<Service name="VirtualEarth">
<ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
</Service>
<MaxConnections>4</MaxConnections>
<Cache/>
</GDAL_WMS>'
gdalio_terra(ve_src, bands = 1:3, band_output_type = "integer")
}
terra::plot(get_im_bing(m))
get_tiles = function(x, ..., plotall = TRUE, native_crs = 27700) {
x_native = sf::st_transform(x, native_crs)
bb_grid_native = sf::st_make_grid(x_native, ...)
bb_grid = sf::st_transform(bb_grid_native, sf::st_crs(x))
bb = sf::st_bbox(bb_grid)
mapview::mapview(bb_grid)
i = 1 # for testing
# browser() # for debugging
im = get_im_bing(bb_grid[i])
if(plotall) terra::plotRGB(im)
i = 2 # for testing
for(i in 1:length(bb_grid)) {
im_n = get_im_bing(bb_grid[i])
r = terra::res(im)
# im_grid = terra::rast(xmin = bb[1], xmax = bb[3], ymin = bb[2], ymax = bb[4], resolution = r) # fail
im_grid = terra::rast(extent = terra::ext(im_n), resolution = terra::res(im))
if(!all(r == terra::res(im_n))) {
im_n = terra::resample(im_n, im_grid)
# terra::res(im_n) = terra::res(im)
}
im = terra::merge(im, im_n)
if(plotall) terra::plotRGB(im)
}
im
}
# im_all = get_tiles(x = m, n = c(6, 3))
im_all = get_tiles(x = m, cellsize = 100)
terra::writeRaster(im_all, "/tmp/im_all.tif")I don't get why you get one overall image, and then do the same by visiting tiles ... get_img_bing() should be enough, also you are using an sf object for an extent but you aren't using the projection of the sf object, or allowing user to set the dimensions
so ...
library(gdalio)
get_im_bing_from_sf = function(x, dimension = c(512, 512), resample = "bilinear") {
bb = sf::st_bbox(x)
ex = as.numeric(bb)[c(1, 3, 2, 4)]
g = gdalio_set_default_grid(list(extent = ex, dimension = dimension, projection = sf::st_crs(x)$wkt))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
ve_src <- '<GDAL_WMS>
<Service name="VirtualEarth">
<ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
</Service>
<MaxConnections>4</MaxConnections>
<Cache/>
</GDAL_WMS>'
gdalio_terra(ve_src, bands = 1:3, band_output_type = "integer", resample = resample)
}
## any sf object will do
ex = c(-2.518027, -2.431069, 52.158692, 52.188771)
m1 = matrix(ex, ncol = 2)
m2 = rbind(m1, m1[1, ])
m = sf::st_polygon(list(m2))
m = sf::st_sfc(m, crs = 4326)
terr <- get_im_bing_from_sf(m)
terra::plotRGB(terr)
the final thing is setting dimension to get the resolution/dimension you want and possibly tweaking the resample. You want dimension to match the aspect ratio of the extent and I'll follow up with that if it's of interest. More stuff might go the gdalio issues, rather than here. Thanks!
no need to conjure up your own crs, just use whatever sf object you have
library(gdalio)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
nc <- read_sf(system.file("shape/nc.shp", package="sf", mustWork = TRUE))
get_im_bing_from_sf = function(x, dimension = c(512, 512), resample = "bilinear") {
bb = sf::st_bbox(x)
ex = as.numeric(bb)[c(1, 3, 2, 4)]
g = gdalio_set_default_grid(list(extent = ex, dimension = dimension, projection = sf::st_crs(x)$wkt))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
ve_src <- '<GDAL_WMS>
<Service name="VirtualEarth">
<ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
</Service>
<MaxConnections>4</MaxConnections>
<Cache/>
</GDAL_WMS>'
gdalio_terra(ve_src, bands = 1:3, band_output_type = "integer", resample = resample)
}
nc <- sf::st_transform(nc, "EPSG:2264")
terr <- get_im_bing_from_sf(nc, dimension = c(1024, 512))
terra::plotRGB(terr)
Created on 2021-10-16 by the reprex package (v2.0.0)
hi to your Da, and congrats on family! ;)
Don't fully understand the code above, I ended up using this:
remotes::install_github("hypertidy/gdalio")
# ex <- c(-2.7171221, -2.7146559, 52.0541044, 52.0549124)
# # get bb interactively
# m = mapedit::editMap()
# mapview::mapview(m)
# bb = sf::st_bbox(m)
# ex = as.numeric(bb)[c(1, 3, 2, 4)]
# dput(ex)
ex = c(-2.518027, -2.431069, 52.158692, 52.188771)
m1 = matrix(ex, ncol = 2)
m2 = rbind(m1, m1[1, ])
m = sf::st_polygon(list(m2))
m = sf::st_sfc(m, crs = 4326)
m_grid = sf::st_make_grid(x = m, n = c(3, 2))
#ex <- c(t(bb))
library(gdalio)
## set target grid (can be anything, any projection), it's just 6 numbers and a string
## can also use stars/terra/raster/grd objects
g <- gdalio_set_default_grid(list(extent = ex, dimension = c(512, 256), projection = "OGC:CRS84"))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
ve_src <- '<GDAL_WMS>
<Service name="VirtualEarth">
<ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
</Service>
<MaxConnections>4</MaxConnections>
<Cache/>
</GDAL_WMS>'
## get image (my _terra _stars _raster wrappers need updating hence band_output_type for now)
im <- gdalio_terra(ve_src, bands = 1:3, band_output_type = "integer")
class(im)
terra::plotRGB(im)
terra::writeRaster(im, filename = "/tmp/cathedral.tif")
get_im_bing = function(x) {
bb = sf::st_bbox(x)
ex = as.numeric(bb)[c(1, 3, 2, 4)]
g = gdalio_set_default_grid(list(extent = ex, dimension = c(512, 256), projection = "OGC:CRS84"))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
ve_src <- '<GDAL_WMS>
<Service name="VirtualEarth">
<ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
</Service>
<MaxConnections>4</MaxConnections>
<Cache/>
</GDAL_WMS>'
gdalio_terra(ve_src, bands = 1:3, band_output_type = "integer")
}
terra::plot(get_im_bing(m))
get_tiles = function(x, ..., plotall = TRUE, native_crs = 27700) {
x_native = sf::st_transform(x, native_crs)
bb_grid_native = sf::st_make_grid(x_native, ...)
bb_grid = sf::st_transform(bb_grid_native, sf::st_crs(x))
bb = sf::st_bbox(bb_grid)
mapview::mapview(bb_grid)
i = 1 # for testing
# browser() # for debugging
im = get_im_bing(bb_grid[i])
# if(plotall) terra::plotRGB(im)
# i = 2 # for testing
# for(i in 1:length(bb_grid)) {
# im_n = get_im_bing(bb_grid[i])
# r = terra::res(im)
# # im_grid = terra::rast(xmin = bb[1], xmax = bb[3], ymin = bb[2], ymax = bb[4], resolution = r) # fail
# im_grid = terra::rast(extent = terra::ext(im_n), resolution = terra::res(im))
# if(!all(r == terra::res(im_n))) {
# im_n = terra::resample(im_n, im_grid)
# # terra::res(im_n) = terra::res(im)
# }
# im = terra::merge(im, im_n)
# f = list.files(tempdir(), pattern = ".tif")
# fs::file_info(f)
# message(f)
# message(file.size(f))
# file.remove(f)
# if(plotall) terra::plotRGB(im)
# }
# im = lapply(bb_grid[1:10], function(x) {
im = lapply(bb_grid, function(x) {
im_n = get_im_bing(x)
im_grid = terra::rast(extent = terra::ext(im_n), resolution = terra::res(im))
if(!all(terra::res(im) == terra::res(im_n))) {
im_n = terra::resample(im_n, im_grid)
# terra::res(im_n) = terra::res(im)
}
im_n
}
)
im
}
# im_all = get_tiles(x = m, n = c(6, 3))
im_all = get_tiles(x = m, cellsize = 100)
saveRDS(im_all, "/tmp/im_all.Rds")
terra::plotRGB(im_all[[1]])
terra::plotRGB(im_all[[2]])
# im_all_flat = terra::merge(im_all[[1]], im_all[[2]], im_all[[3]])
im_all_flat = do.call(terra::merge, im_all)
dir.create("tifs")
i = 1
for(i in 2:length(im_all)) {
i_c = formatC(i, digits = 3, flag = 0)
message(i_c)
f = file.path("tifs", paste0(i_c, ".tif"))
terra::writeRaster(im_all[[i]], f)
message(file.size(f))
}
system("gdalbuildvrt mosaic.vrt tifs/*.tif")
system('gdal_translate -of GTiff -co "TILED=YES" mosaic.vrt mosaic.tif')
system('curl --upload-file mosaic.tif https://transfer.sh/ ')
system('gdal_translate -of GTiff -co NBITS -co "COMPRESS=JPEG" -co "TILED=YES" mosaic.vrt mosaic2.tif')
zip(zipfile = "mosaic.zip", files = "mosaic.tif")
piggyback::pb_upload("mosaic.tif")
terra::plotRGB(im_all_flat)
terra::writeRaster(im_all, "/tmp/im_all.tif")
I don't know what else I can say ... set the grid you want, run gdalio_terra - all the rest is redoing the same task a more complicated way
But does the service not give you lower res data if you set a large grid? I should try it. V. inefficient code 🤯 and doh! if not!
no, it gives what you ask for 😉 I know I'm butting up against lots of established habit ... but it really is much simpler this way
... there are limits, and it doesn't yet support out of memory for really large jobs, but a decent image is easy