How to access WMS Servers in R Programming Language

How to access WMS Servers in R Programming Language

I couldn’t find a good package to read WMS data into a SpatialRaster or for downloading. However I discovered that you can use sf’s built-in version of gdal_util query the WMS server’s information then use the sf built in gdal_translate to download the WMS image into a temporary file then load it into memory.

view source
print?
01library(sf)
02library(tidyverse)
03library(terra)
04rm(list=ls())
05 
06# obtain layer info from WMS server using sf built-in version of gdal_info
07wms_url <- 'https://orthos.its.ny.gov/ArcGIS/services/wms/Latest/MapServer/WMSServer?'
08ginfo <- sf::gdal_utils('info', str_c('WMS:',wms_url), quiet=T)
09 
10# extract layers and layer urls from returned layer info,
11# create a data table from this data
12ldesc <- ginfo %>% str_match_all('SUBDATASET_(\\d)_DESC=(.*?)\n')
13ldsc <- ldesc[[1]][,3]
14lurl <- ginfo %>% str_match_all('SUBDATASET_(\\d)_NAME=(.*?)\n')
15lurl <- lurl[[1]][,3]
16wms_layers <- cbind(ldesc, lurl)
17rm(ldesc, lurl)
18 
19# obtain from census bureau a shapefile of the Town of Cazenovia,
20# then find the bounding box around it
21bbx<- tigris::county_subdivisions('ny') %>%
22  filter(NAME == 'Cazenovia') %>%
23  st_transform(3857) %>%
24  st_bbox()
25 
26# Revise the WMS URL based on the above bounding box
27# Web Mercator (3857) is best to use with WMS Servers as
28# that is the native projection. As I wanted multiple layers
29# I also revised the URL string to include all layers desired
30# Seperated by a comma (explore the wms_layer table for details)
31url <- wms_layers[4,2] %>%
32  str_replace('LAYERS=(.*?)&', 'LAYERS=3,2,1,0&') %>%
33  str_replace('SRS=(.*?)&', 'SRS=EPSG:3857&') %>%
34  str_replace('BBOX=(.*?)$', str_c('BBOX=',paste(bbx, collapse=',')) )
35 
36# Use sf built-in gdal_utils to download the image of Cazenovia based
37# on the URL, with an output size listed below. You can also download
38# based on resolution by using the tr option that is commented out below.
39# then load the temporary file into rast as a Spatial Raster for further
40# processing. In addition, adding -co COMPRESS=JPEG greatly reduces the size
41# of aerial photography with minimal impacts on quality or loading speed.
42 
43t <- tempfile()
44reso <-  c('-outsize','1920','1080','-co','COMPRESS=JPEG')
45#reso <- c('-tr', '1','1','-co','COMPRESS=JPEG') # "meters" in the 3857 projection
46sf::gdal_utils('translate', url, t, reso)
47r <- rast(t)
48 
49# display the spatial raster or whatever you would like to do with it
50# the tiff is stored in the temporary location in the t variable
51plot(r)

Leave a Reply

Your email address will not be published. Required fields are marked *