Skip to content

Commit

Permalink
master #4 - correct indent., add findP4s to namespace, fix urn regexp
Browse files Browse the repository at this point in the history
  • Loading branch information
eblondel committed Nov 21, 2013
1 parent 49066df commit d5059bd
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 65 deletions.
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ import(rgeos)
import(dplR)

export(
findP4s,
readWFS,
exportFeatures,
getIntersection
Expand Down
130 changes: 65 additions & 65 deletions R/FigisGeoUtils.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down

0 comments on commit d5059bd

Please sign in to comment.