diff --git a/NAMESPACE b/NAMESPACE index 312e2ad..349b44d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,6 +5,7 @@ import(rgeos) import(dplR) export( + findP4s, readWFS, exportFeatures, getIntersection diff --git a/R/FigisGeoUtils.R b/R/FigisGeoUtils.R index 2a2e259..bc7fdb1 100644 --- a/R/FigisGeoUtils.R +++ b/R/FigisGeoUtils.R @@ -18,38 +18,38 @@ # We actually except lat/long coordinates systems due to confusion in output of different WFS vendors # findP4s <- function(srsName, morphToESRI=FALSE) { - - if (missing(srsName)) { - stop("please provide a spatial reference system name") - } - #we remove the latlong proj for compatibility with sp - proj.lst <- as.character(projInfo("proj")$name) - proj.lst <- proj.lst[proj.lst != "latlong" & proj.lst != "latlon"] - #build combinations of know proj and datum - proj.datum.grd <- expand.grid(proj=proj.lst, datum=as.character(projInfo("datum")$name), stringsAsFactors=FALSE) - #function to ask WKT representation - getShowWkt <- function(x) { - res <- try(showWKT(paste("+proj=", x[1], " +datum=", x[2], sep=""), morphToESRI=morphToESRI), silent=TRUE) - if (class(res) == "try-error") { - return(NA) - } else { - return(res) - } - } - #ask WKT for all projection - GCS <- apply(proj.datum.grd, 1, FUN=getShowWkt) - - GCS.df <- data.frame(proj=proj.datum.grd$proj, datum=proj.datum.grd$datum, WKT=GCS, stringsAsFactors=FALSE) - #keep only valids - GCS.df <- GCS.df[! is.na(GCS.df$WKT),] - #the pattern to find - pattern <- paste("GEOGCS[\"", srsName, "\"", sep="") - #search for pattern - GCS.df <- GCS.df[substr(tolower(GCS.df$WKT), 1, nchar(pattern)) == tolower(pattern),] - #keep only first SRS in case of identical WKT representation - GCS.df <- GCS.df[!duplicated(GCS.df$WKT),] - #return the proj4 definition - return(paste("+proj=", GCS.df$proj, " +datum=", GCS.df$datum, sep="")) + + if (missing(srsName)) { + stop("please provide a spatial reference system name") + } + #we remove the latlong proj for compatibility with sp + proj.lst <- as.character(projInfo("proj")$name) + proj.lst <- proj.lst[proj.lst != "latlong" & proj.lst != "latlon"] + #build combinations of know proj and datum + proj.datum.grd <- expand.grid(proj=proj.lst, datum=as.character(projInfo("datum")$name), stringsAsFactors=FALSE) + #function to ask WKT representation + getShowWkt <- function(x) { + res <- try(showWKT(paste("+proj=", x[1], " +datum=", x[2], sep=""), morphToESRI=morphToESRI), silent=TRUE) + if (class(res) == "try-error") { + return(NA) + } else { + return(res) + } + } + #ask WKT for all projection + GCS <- apply(proj.datum.grd, 1, FUN=getShowWkt) + + GCS.df <- data.frame(proj=proj.datum.grd$proj, datum=proj.datum.grd$datum, WKT=GCS, stringsAsFactors=FALSE) + #keep only valids + GCS.df <- GCS.df[! is.na(GCS.df$WKT),] + #the pattern to find + pattern <- paste("GEOGCS[\"", srsName, "\"", sep="") + #search for pattern + GCS.df <- GCS.df[substr(tolower(GCS.df$WKT), 1, nchar(pattern)) == tolower(pattern),] + #keep only first SRS in case of identical WKT representation + GCS.df <- GCS.df[!duplicated(GCS.df$WKT),] + #return the proj4 definition + return(paste("+proj=", GCS.df$proj, " +datum=", GCS.df$datum, sep="")) } # Read WFS & returns a sp object @@ -68,7 +68,7 @@ readWFS <- function(url, outputFormat = "GML", p4s = NULL, gmlIdAttributeName="g if(outputFormat != "GML") { wfsRequest <- paste(url, "&outputFormat=", outputFormat, sep="") } - + if(outputFormat == "GML") { # download the data @@ -77,43 +77,43 @@ readWFS <- function(url, outputFormat = "GML", p4s = NULL, gmlIdAttributeName="g download.file(wfsRequest, destfile, mode="wb") layername <- ogrListLayers(destfile) if (length(layername) != 1) { - stop("Mistake with layers in the input dataset") + stop("Mistake with layers in the input dataset") } - + # get the Spatial Reference System (SRS) srs <- NA xmlfile<-xmlTreeParse(destfile, useInternalNodes = TRUE) srsName <- getNodeSet(xmlfile, "(//gml:featureMember//@srsName)[1]") - - 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: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 - srs <- findP4s(srsName, morphToESRI=FALSE) - if (length(srs) == 0) { - #if not found search with consideration of the ESRI representation - srs <- findP4s(srsName, morphToESRI=TRUE) - } - if (! length(srs) == 1) { - srs <- NA - } - } - } - } + + 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 + srs <- findP4s(srsName, morphToESRI=FALSE) + if (length(srs) == 0) { + #if not found search with consideration of the ESRI representation + srs <- findP4s(srsName, morphToESRI=TRUE) + } + if (! length(srs) == 1) { + srs <- NA + } + } + } + } if(is.na(srs)){ warning("Unable to convert GML srsName to a CRS object. CRS will be set to NA", call. = T)