R/readWFS.R

In openfigis/RFigisGeo: Set of geoprocessing facilities used in FIGIS

#' @name readWFS
#' @title readWFS
#' @description Read WFS & returns a sp object
#'
#' @param url the baseURL of the WFS GetFeature request
#' @param outputFormat the output format for the WFS GetFeature request, by default "GML"
#' @param p4s an optional proj4string, by default NULL (an attempt will be performed to get the projection from the data)
#' @param gmlIdAttributeName specific to GML, the name of the ID attribute, by default "gml_id"
#' @param target.dir a target directory where temporary GML files will be downloaded
#' @param verbose if log has to printed in the R console. TRUE by default
#' @return an object of class "Spatial"
#' 
#' @author Emmanuel Blondel \email{emmanuel.blondel1@@gmail.com}
#'         Norbert Billet \email{norbert.billet@@ird.fr} 
#
readWFS <- function(url, outputFormat = "GML", p4s = NULL,
                    gmlIdAttributeName="gml_id", target.dir = NULL,
                    verbose = TRUE){
  
  logger.warn(paste("The 'readWFS' function is deprecated. Please use 'ows4R' package (available on CRAN) with the 'WFSClient'.",
                   "Examples are available at https://github.com/eblondel/ows4R/wiki#43-r-interface-to-web-feature-service-wfs"))

  features <- NULL
  
  #request
  wfsRequest <- url
  if(outputFormat != "GML") {
    wfsRequest <- paste(url, "&outputFormat=", outputFormat, sep="")
  }
  
  if(verbose) logger.info(sprintf("Reading WFS request: %s \n",wfsRequest))
  if(outputFormat == "GML") {
    # download the data
    
    #escaping characters in cql_filter
    percentCodes <- c("%20", "%27")
    names(percentCodes) <- c("\x20", "\x27")
    if(grepl("cql_filter", wfsRequest)){
      arr <- unlist(strsplit(wfsRequest, "&cql_filter="))
      wfsBaseRequest <- arr[1]
      cql_filter <- arr[2]
      otherParams <- NULL
      arr_suffix <- unlist(strsplit(cql_filter, "&"))
      if(length(arr_suffix)>1){
        cql_filter <- arr_suffix[1]
        otherParams <- arr_suffix[2:length(arr_suffix)]
      }
      cql_filter <- curlPercentEncode(x = curlEscape(cql_filter), codes = percentCodes)
      wfsRequest <- paste0(wfsBaseRequest, "&cql_filter=", cql_filter)
      if(!is.null(otherParams)){
        wfsRequest <- paste0(wfsRequest, "&", paste(otherParams, collapse="&"))
      }
    }
    
    content <- getURL(wfsRequest, .encoding = "UTF-8")
    xmlfile <- xmlTreeParse(content, useInternalNodes = TRUE)
    
    if(!all(class(xmlfile) == c("XMLInternalDocument","XMLAbstractDocument"))){
      if(verbose) logger.error("Fetch data doesn't seem in XML format \n")
      return(NULL)
    }
    
    #write the file to disk
    tempdir <- NULL
    if(missing(target.dir) || is.null(target.dir)){
      tempdir <- tempdir()
    }else{
      if(file.exists(target.dir)){ 
        tempdir <- target.dir
      }else{
        if(verbose) logger.error(sprintf("Target directory %s doesn't exist \n",target.dir))
        return(NULL)
      }
    }
    
    if(verbose){
      logger.info(sprintf("Writing temporary GML file to '%s' \n",tempdir))
    }
    
    tempf = tempfile(tmpdir = tempdir) 
    destfile = paste(tempf,".gml",sep='')
    saveXML(xmlfile, destfile)
    
    if(!file.exists(destfile)){
      if(verbose) logger.error(sprintf("GML file %s cannot be found \n", destfile))
      return(NULL)
    }
    
    deleteGML <- function(){
      if(verbose) logger.info(sprintf("Deleting temporary GML file '%s' \n",destfile))
      unlink(destfile)
      unlink(paste0(tempf,".gfs"))
    }
    
    #download.file(wfsRequest, destfile, mode="wb")
    if(verbose) logger.info(sprintf("Fetching layer info from GML file '%s' using GDAL \n",destfile))
    layername <- tryCatch(
      ogrListLayers(destfile),
      error = function(error) {
        if(verbose){
          logger.error(error)
          deleteGML()
        }
      })
    
    
    if(is.null(layername) || length(layername) == 0) {
      if(verbose) logger.error("Unknown or Empty GIS web-resource \n")
      return(NULL)
    }else{
      if(verbose) logger.info(sprintf("Layer identified by GDAL as '%s'", layername))
    }
    
    #check if we have geometry
    propertyNames <- NULL
    properties <- unlist(strsplit(wfsRequest,"&propertyName="))
    if(length(properties) > 1){
      propertyNames <- unlist(strsplit(properties[length(properties)], "&"))[1]
      propertyNames <- unlist(strsplit(propertyNames, ","))
    }
    if(!is.null(propertyNames)){
      hasGeometry <- any(c("GEOMETRY", "geometry", "THE_GEOM","the_geom") %in% propertyNames)
    }else{
      hasGeometry = ((length(getNodeSet(xmlfile, "//gml:featureMember//gml:coordinates")) > 0)
                     || (length(getNodeSet(xmlfile, "//gml:featureMember//gml:pos")) > 0)
                     || (length(getNodeSet(xmlfile, "//gml:featureMember//gml:posList")) > 0))
    }
    
    if(hasGeometry){
      
      # get the Spatial Reference System (SRS)
      if(verbose) logger.info("Inheriting Spatial Reference System...\n")
      srs <- NA
      #xmlfile<-xmlTreeParse(destfile, useInternalNodes = TRUE)
      fmXML <- getNodeSet(xmlfile, "//gml:featureMember")
      if(length(fmXML) > 0)
        fmXML <- fmXML[[1]]
      srsName <- getNodeSet(xmlDoc(fmXML), "//@srsName")
      if (length(srsName) >= 1) {
        srsName <- as.character(srsName[[1]])
      }	
      #srsName patterns matching		
      srsPattern = "http://www.opengis.net/gml/srs/epsg.xml#" #match case 1
      if(attr(regexpr(srsPattern, srsName, ignore.case = T),"match.length") > 0){
        epsg <- unlist(strsplit(srsName, srsPattern))[2]
        srs <- paste("+init=epsg:", epsg, sep="")
      }else{
        srsPattern = "urn:(x-)?ogc:def:crs:EPSG" #match case 2
        if(attr(regexpr(srsPattern, srsName, ignore.case = T),"match.length") > 0){
          srsStr <- unlist(strsplit(srsName, ":"))
          epsg <- srsStr[length(srsStr)]
          srs <- paste("+init=epsg:", epsg, sep="")
        }else{
          #search if srsName is a WKT PROJ name (PROJCS or GEOGCS)
          #if yes set srs with the corresponding proj4string
          #first search without any consideration of the ESRI representation
          if(verbose) logger.info("Fetching SRS definition\n")
          srs <- findP4s(srsName, morphToESRI=FALSE)
          if (is.na(srs)){
            #if not found search with consideration of the ESRI representation
            srs <- findP4s(srsName, morphToESRI=TRUE)
          }
          if (! is.na(srs) && ! length(srs) == 1) {
            srs <- NA
          }
        }
      }
      
      if(is.na(srs)){
        if(verbose) logger.warn("Unable to convert GML srsName to a CRS object. CRS will be set to NA \n")
      }else{
        if(verbose) logger.info(sprintf("SRS definition = '%s' \n",srs))
      }
      
      if (missing(p4s)) p4s <- srs
      if(verbose) logger.info(sprintf("Reading temporary GML file '%s' with GDAL \n",destfile))
      features = tryCatch(
        readOGR(destfile, layername, p4s = srs, disambiguateFIDs = TRUE, verbose = verbose),
        error = function(err){
          if(verbose){
            logger.error(err)
            deleteGML()
          }
        })
      if(!is.null(features)){
        if(regexpr("SpatialPoints", class(features)) == -1){
          logger.info(sprintf("Add feature identifiers for object of class '%s' \n", class(features)))
          
          if(!(gmlIdAttributeName %in% colnames(slot(features, "data")))){
            slot(features, "data") <- cbind(gml_id = paste(layername,1:nrow(slot(features,"data")),sep="."), slot(features,"data"))
          }
          features <- spChFIDs(features, as.character(features@data[,gmlIdAttributeName]))
        }
      }else{
        if(verbose) logger.warn(sprintf("Features returned by GDAL are null \n",destfile))	
      }
      
    }else{
      featureMembers <- getNodeSet(xmlfile, "//gml:featureMember")
      if(length(featureMembers)>0){
        membersContent <- sapply(featureMembers, function(x) xmlChildren(x))
        fid <- sapply(membersContent, function(x) xmlAttrs(x))
        membersAttributes <- xmlToDataFrame(
          nodes = getNodeSet(xmlfile, "//gml:featureMember/*[@*]"),
          stringsAsFactors = FALSE
        )
        features <- cbind(fid, membersAttributes, stringsAsFactors = FALSE)
      }
    }
    
    #unlink temporary files
    deleteGML()
    
  }else{
    stop("Unsupported WFS format")
  }
  
  if(verbose){
    if(!is.null(features)){
      logger.info("WFS features successfully fetched! \n")
      logger.info(sprintf("Number of features = %s \n", nrow(features)))
    }else{
      logger.warn("Empty WFS feature collection \n")
    }
  }
  return(features)
}