Extracting data from WMS using R

wms_data.md

WMS in R from https://gis.stackexchange.com/questions/416153/extracting-data-from-wms-using-r-is-it-possible

and https://twitter.com/TimSalabim3/status/1458856094281916419?s=20

##I chose one from these options
vapour::vapour_sds_names("https://www.statistik.at/gs-open/ATLAS_UNFALL_OPEN/wms?service=WMS&version=1.1.0&request=GetCapabilities")$subdataset
#> WMS ...
#> ...
#> ...
#> 

## to obtain this dsn use gdalinfo or sds_names as above, might need /vsicurl ahead of the sds/gdalinfo query to avoid driver conflation
dsn <- "WMS:http://www.statistik.at/gs-open/ATLAS_UNFALL_OPEN/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=unfall_koord&SRS=EPSG:3857&BBOX=1061505.0695717288,5840291.6646587765,1908795.0257996104,6274985.801868425"

## get info on the source
(ri <- vapour::vapour_raster_info(dsn))
#> $geotransform
#> [1]  1.061505e+06  7.891003e-04  0.000000e+00  6.274986e+06  0.000000e+00
#> [6] -7.891003e-04
#> 
#> $dimXY
#> [1] 1073741824  550873136
#> 
#> $minmax
#> [1] NA NA
#> 
#> $tilesXY
#> [1] 1024 1024
#> 
#> $projection
#> [1] "PROJCS[\"WGS 84 / Pseudo-Mercator\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],EXTENSION[\"PROJ4\",\"+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\"],AUTHORITY[\"EPSG\",\"3857\"]]"
#> 
#> $bands
#> [1] 3
#> 
#> $projstring
#> [1] "+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"
#> 
#> $nodata_value
#> [1] -1e+10
#> 
#> $overviews
#>  [1] 536870912 275436568 268435456 137718284 134217728  68859142  67108864
#>  [8]  34429571  33554432  17214786  16777216   8607393   8388608   4303696
#> [15]   4194304   2151848   2097152   1075924   1048576    537962    524288
#> [22]    268981    262144    134491    131072     67245     65536     33623
#> [29]     32768     16811     16384      8406      8192      4203      4096
#> [36]      2101      2048      1051      1024       525
#> 
#> $filelist
#> [1] NA
#> 
#> $datatype
#> [1] "Byte"
#> 
#> $extent
#> [1] 1061505 1908795 5840292 6274986
# remotes::install_github("hypertidy/gdalio")
library(gdalio)
## set up grid based on source info, but choose a low-res overview for the dimension
gdalio_set_default_grid(list(extent = ri$extent, dimension = tail(ri$overviews, 4)[1:2], 
                             projection = ri$projection))
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))

terra::plotRGB(gdalio_terra(dsn, bands = 1:3))

^{Created on 2021-11-12 by the reprex package (v2.0.1)}

Hi, I know this wasn't my original question at all, but now just curious. To increase the resolution, you just increase the dimension variable right? Is there an upper bound to that, 536870912 275436568?

How would you interpret the bands? I know there different spatraster layers. The xml files indicates that the layer has variables on for car, bike, pedestrian, other accidents, and date and location variables, so wondering what is actually being graphed. Thanks!

there's no upper bound but no computer known will handle the result 😀, the native resolution upper bound is in the '$overviews', those are pairs of dimensions cascaded down from '$dimXY', they are the dimensions at each xoom level

bands are just red green blue as per '$bands', but different layers are in the subdataset list, you have to query and read from those separately - see 'LAYERS=unfall_koord', the full list of available sds comes from 'vapour_sds_names', I just commented that out

note that you don't have to use dimensions or projection or extent from the source, you can use any set of these for a map projection

but, sensibly limited to an in-memory request - there is no handling of very large reads here, you would use this to make a number of smaller requests for that