How to save WMS images to raster with ows4R so that they can be plotted with other packages such as tmap?

[Florent-Demoraes]

Hello,

I was wondering wether they would be a way of downloading images from a WMS server so as to create a local raster file and to plot it with other packages such as tmap for instance.

Here is what I am doing til now, but it is quite complicated and it does not rely on ows4R package.

The code below is derived from that published here : https://gist.github.com/obrl-soil/73916c6fe223d510293cb1a2bbf2879a

library(sf)
library(raster)
library(httr)
library(tmap)
library(RStoolbox)

web service - QLD Imagery Whole of State Public basemap (no WCS unforts, hence all this)

qldimg_wms <- 'https://spatial-img.information.qld.gov.au/arcgis/services/Basemaps/LatestSatelliteWOS_AllUsers/ImageServer/WMSServer?'

Cubesat 2.4m res, from Q3 2017 @ Oct 2019

going to grab a bit from over Brisbane CBD for testing.

bb <- structure(c(153.00795, -27.49034, 153.04041, -27.45917),
names = c("xmin", "ymin", "xmax", "ymax"),
class = "bbox", crs = sf::st_crs(4326))

bbstr <- paste0(as.vector(bb), collapse = ',')

bbm <- st_bbox(st_transform(st_as_sfc(bb), 28356)) # local UTM grid

about as close to src alignment as u need for a quick grab:

w <- plyr::round_any(abs(bbm[[3]] - bbm[[1]]), 2.4) / 2.4
h <- plyr::round_any(abs(bbm[[2]] - bbm[[4]]), 2.4) / 2.4

construct WMS URL

img_url <-
paste0(qldimg_wms, 'service=WMS', '&version=1.3.0', '&request=GetMap',
'&layers=LatestSatelliteWOS_AllUsers', '&styles=', '&crs=CRS%3A84',
'&bbox=', bbstr, '&width=', w, '&height=', h, '&format=image%2Fpng')

download and save locally

httr::GET(url = img_url
, write_disk(file.path(getwd(), 'BNE.png'), overwrite = TRUE)
)

import PNG

img <- png::readPNG(file.path(getwd(), 'BNE.png'))

convert each band to a raster

rasR <- raster(img[,,1], xmn = bb[1], xmx = bb[3], ymn = bb[2], ymx = bb[4])
rasG <- raster(img[,,2], xmn = bb[1], xmx = bb[3], ymn = bb[2], ymx = bb[4])
rasB <- raster(img[,,3], xmn = bb[1], xmx = bb[3], ymn = bb[2], ymx = bb[4])

create a stack

ras <- stack(rasR, rasG, rasB)

georeference image

crs(ras) <- CRS('+init=EPSG:4326')

write image

writeRaster(ras, file.path(getwd(), 'BNE.tif'))

plotting options:

plotRGB in base

plotRGB(ras, scale = 1, maxpixels = ncell(ras), interpolate = TRUE)

tm_rgb in tmap (bit slow but has interp = T by default)

tm_shape(ras) +
tm_rgb(max.value = 1)

ggRGB in RStoolbox (ggplot-friendly wrapper for the above)

ggRGB(ras, r = 1, g = 2, b = 3, maxpixels = ncell(ras))

[eblondel]

@Florent-Demoraes thanks for this, what you shared is actually really useful for investigating how to plug the getMap operation. As soon as I have some time, I will start investigating it

[paul-carteron]

Hi,
If that can help, to implement getMap operation in happign::get_wms_raster() I used gdal warp operation from sf::gdal_utils().

It's quite straightforward :

• Building url for GDAL with WMS driver
• Use gdal_utils("warp") to a tempfile
• Import in R with terra::rast()