Skip to content

Instantly share code, notes, and snippets.

@jlehtoma
Created November 13, 2011 20:01

Revisions

  1. jlehtoma created this gist Nov 13, 2011.
    176 changes: 176 additions & 0 deletions wms_proto.R
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,176 @@
    # Author: Joona Lehtomäki <joona.lehtomaki@gmail.com>
    # Updated: 13.11.2011
    # Version: 0.0.1

    if (!require("rgdal")) {
    install.packages("rgdal")
    }

    if (!require("raster")) {
    install.packages("raster")
    }

    if (!require("XML")) {
    install.packages("XML")
    }

    if (!require("sorvi")) {
    install.packages("sorvi", repos="http://R-Forge.R-project.org",
    type = "source", dependencies = TRUE)
    }

    # URLs to OIVA site (http://wwwp2.ymparisto.fi/scripts/oiva.asp) WMS servers.
    # Servers registered here include two data sets so far:
    # 1. Maanpeite (Corine 2000 / 2006) - land use categorization
    # 2. Suojelualueet (suojelualueet + Natura2000) - protected areas
    OIVA.URLS <- c(Corine='http://paikkatieto.ymparisto.fi/ArcGIS/services/INSPIRE/SYKE_Maanpeite/MapServer/WMSServer?',
    Suojelu='http://paikkatieto.ymparisto.fi/ArcGIS/services/INSPIRE/SYKE_SuojellutAlueet/MapServer/WMSServer?')

    #' Create a WMS service description that is used to load data from WMS server.
    #'
    #' WMS service description file is a XML file that describes required and
    #' optional information on how to retrieve an exisiting WMS raster over the
    #' web. The extent of the raster tile from the data source is defined by the
    #' extent of a SpatialPolygonsDataFrame object (no other ways of
    #' providing extent are implemented yet). Raster resolution (pixel size is
    #' also provided as a parameter.
    #'
    #' @param data.source string describing the name of the data source
    #' @param layer string name of the layer to be fetched from the data source
    #' @param extent SpatialPolygonsDataFrame object to be used to define the extent
    #' @param resolution integer value of the resolution (CRS dependent)
    #'
    #' @return A character string of XML
    #'
    #' @note function does not check whether data source or the layer actually exist
    #'
    #' @references
    #' \url{http://www.gdal.org/frmt_wms.html}
    #' @author Joona Lehtomäki \email{joona.lehtomaki@@gmail.org}

    WMSserviceDesc <- function(data.source, layer, extent, resolution) {

    # Check if data.source is included in OIVA.URLS
    if (is.na(OIVA.URLS[data.source])) {
    stop(paste('Data source', data.source, 'not found'))
    }

    # Extent is defined by the bounding box of the SpatialPolygonsObject provided
    # as extent parameter.
    # TODO: implement other ways of providing the extent
    if (class(extent) == 'SpatialPolygonsDataFrame') {
    bbox.extent <- bbox(extent)
    # Number of columns and rows (i.e. resolution) is defined by the real
    # width and height of the raster divided by the resolution parameter
    # (all this depends on the CRS, not very tested)
    ncols <- round((bbox.extent[1, 2] - bbox.extent[1, 1]) / resolution, 0)
    nrows <- round((bbox.extent[2, 2] - bbox.extent[2, 1]) / resolution, 0)
    # Set the extent corners
    ulx <- bbox.extent[1, 1]
    uly <- bbox.extent[2, 2]
    lrx <- bbox.extent[1, 2]
    lry <- bbox.extent[2, 1]
    } else {
    stop('Function only supports SpatialPolygonDataFrames')
    }

    # Create the XML structure based on the GDAL specification at
    # http://www.gdal.org/frmt_wms.html

    # Root level
    root <- newXMLNode('GDAL_WMS')
    # Service level
    service <- newXMLNode('Service', attrs=c(name='WMS'), parent=root)
    version <- newXMLNode('Version', text='1.1.1', parent=service)
    serverUrl <- newXMLNode('ServerUrl', text=OIVA.URLS[data.source], parent=service)
    # TODO: CRS should not be hard coded
    srs <- newXMLNode('SRS', text='EPSG:3067', parent=service)
    # Not sure if really needed
    imageFormat <- newXMLNode('ImageFormat', text='image/tiff', parent=service)
    layers <- newXMLNode('Layers', text=layer, parent=service)
    # Style is needed even if empty
    style <- newXMLNode('Style', parent=service)

    # DataWindow level
    dataWindow <- newXMLNode('DataWindow', parent=root)
    # Note that the following notation is minX, maxY, maxX, minY
    upperLeftX <- newXMLNode('UpperLeftX', text=ulx, parent=dataWindow)
    upperLeftY <- newXMLNode('UpperLeftY', text=uly, parent=dataWindow)
    lowerRightX <- newXMLNode('LowerRightX', text=lrx, parent=dataWindow)
    lowerRightY <- newXMLNode('LowerRightY', text=lry, parent=dataWindow)
    # TODO: although size is set here, it is not completely clear how the
    # native raster resolution on the WMS server is related to resolution
    # requested
    sizeX <- newXMLNode('SizeX', text=ncols, parent=dataWindow)
    sizeY <- newXMLNode('Sizey', text=nrows, parent=dataWindow)
    yOrigin <- newXMLNode('YOrigin', text='top', parent=dataWindow)

    # Back to the root level
    projection <- newXMLNode('Projection', text='EPSG:3067', parent=root)
    # Optional, this is also the default. Seems to be required in case where the
    # the raster requested is 3-band RGB raster.
    bandsCount <- newXMLNode('BandsCount', text='3', parent=root)
    # Optional, probably not needed here
    cache <- newXMLNode('Cache', parent=root)

    # Save the created XML object, not providing a file path converts the object
    # into string.
    return(saveXML(root))
    }

    # Use soRvi to get some reference data (municipalities)
    data(MML)
    sp <- MML[["1_milj_Shape_etrs_shape"]][["kunta1_p"]]

    # Preprocess MML Shape file
    sp <- preprocess.shape.mml(sp)

    # Select 2 example locations: Helsinki and Tampere

    # Helsinki - Corine 2006 land use classes

    # Extract only the polygon for Helsinki
    sp.helsinki <- sp[which(sp@data$Kunta_ni1 == "Helsinki"),]
    # Create the WMS description XML string
    wms.xmlDesc.helsinki <- WMSserviceDesc('Corine', 'CorineLandCover2006_25m',
    sp.helsinki, 25)

    # Use GDAL to read in the value, readGDAL returns a SpatialObject
    wms.img.helsinki <-readGDAL(wms.xmlDesc.helsinki)
    # Use brick function from package raster to squeeze the 3 RGB band into a
    # RasterBrick object for easy plotting
    b.wms.img.helsinki <- brick(wms.img.helsinki)
    # Use plotRGB from package raster to plot the RGB raster
    plotRGB(b.wms.img.helsinki)
    # Add municipality borders
    plot(sp.helsinki, add=TRUE, lwd=2)

    # Helsinki - Natura 2000 polygons
    wms.xmlDesc.helsinki <- WMSserviceDesc('Suojelu',
    'ProtectedSites.Natura2000Polygons',
    sp.helsinki, 25)

    wms.img.helsinki <- readGDAL(wms.xmlDesc.helsinki)
    b.wms.img.helsinki <- brick(wms.img.helsinki)
    plotRGB(b.wms.img.helsinki)
    plot(sp.helsinki, add=TRUE, lwd=2)

    # Tampere - Corine 2006
    sp.tampere <- sp[which(sp@data$Kunta_ni1 == "Tampere"),]
    wms.xmlDesc.tampere <- WMSserviceDesc('Corine', 'CorineLandCover2006_25m',
    sp.tampere, 25)

    wms.img.tampere <- readGDAL(wms.xmlDesc.tampere)
    b.wms.img.tampere <- brick(wms.img.tampere)
    plotRGB(b.wms.img.tampere)
    plot(sp.tampere, add=TRUE, lwd=2)

    # Tampere - Natura 2000 polygons
    wms.xmlDesc.tampere <- WMSserviceDesc('Suojelu',
    'ProtectedSites.Natura2000Polygons',
    sp.tampere, 25)

    wms.img.tampere <- readGDAL(wms.xmlDesc.tampere)
    b.wms.img.tampere <- brick(wms.img.tampere)
    plotRGB(b.wms.img.tampere)
    plot(sp.tampere, add=TRUE, lwd=2)