How to extract and georeference satellite WMS imagery in R

satimg_WMS_R

# Recipe for pulling imagery crops from the Queensland Government's public imagery basemap
# remember to credit output images appropriately using the metadata info at
# http://qldspatial.information.qld.gov.au/catalogue/custom/search.page?q=%22Queensland+Imagery+Whole+Of+State+Satellite+Public+Basemap+Service%22

library(sf)
library(raster)
library(httr)

# 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 (don't worry about alpha band)
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])

# stack
ras <- stack(rasR, rasG, rasB)

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

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

# plotting options:

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

# ggRGB in RStoolbox (ggplot-friendly wrapper for the above)
ggRGB(ras, r = 1, g = 2, b = 3, maxpixels = ncell(ras)) +
  theme_minimal() +
  theme(axis.title = element_blank())

# tm_rgb in tmap (bit slow but has interp = T by default)
tm_shape(ras) +
  tm_rgb(max.value = 1)

# NB default colours in QGIS will look a little odd until you adjust the minmax
# values on each band to [0, 1]