Nothing
'.prepare.con' <- function(x,...)
{
arglist <- list(...)
if (missing(x))
x <- NA
if ((!missing(x))&&(!is.ursa(x)))
{
arglist$temporal <- x
ind <- .grep("^temporal$",names(arglist))
names(arglist)[ind] <- ""
x <- .raster.skeleton()
x$grid <- session_grid()
if (.is.grid(arglist$temporal))
x$grid <- arglist$temporal
}
isCT <- FALSE
ind <- .grep("^(colortable|category)",names(arglist))
if (length(ind))
{
ct <- arglist[[ind]]
if (is.logical(ct)) {
if (ct)
ct <- ursa_colortable(x)
else {
ct <- character(0)
ursa_colortable(x) <- ct
}
}
if ((is.character(ct))&&(length(ct))) {
isCT <- TRUE
n <- length(ct)
arglist$datatype <- ifelse(n<255,1L,ifelse(n<65535,12L,13L))
arglist$nodata <- n
class(x$value) <- "ursaCategory"
x$colortable <- ct
class(x$colortable) <- "ursaColorTable"
}
}
simple <- attr(x,"copyright")
simple <- ((is.character(simple))&&(nchar(simple)>0))
if ((simple)) {
if ((!FALSE)||("ursa" %in% loadedNamespaces())) {
g1 <- ursa(x,"grid")
ursa(x,"grid") <- regrid(g1,setbound=c(0,0,g1$columns,g1$rows),crs="")
}
}
x$con <- .create.con(x,arglist)
if (!.is.con(x$con)) {
cat("Unable to create connection\n")
return(NULL)
}
ind <- .grep("^(band|layer)*name",names(arglist))
if (length(ind))
x$name <- arglist[[ind[1]]]
if ((0)&&(length(ind)))
cat(sprintf("x$name=%s\n",x$name))
if (any(is.na(x$dim)))
x$dim <- as.integer(with(x$con,c(lines*samples,bands)))
if ((x$dim[2]!=x$con$bands)&&((length(x$value)==1)&&(is.na(x$value))||(x$dim[2]!=dim(x$value)[2])))
x$dim[2] <- x$con$bands
if (length(.grep("(band|layer|names)",names(arglist))))
{
if (x$dim[2]!=x$con$bands)
x$dim[2] <- x$con$bands
}
else
{
if (x$con$bands!=x$dim[2])
x$con$bands <- x$dim[2]
}
if ((TRUE)&&(!is.na(x$name[1]))&&(!is.na(x$con$posZ[1]))) ## FALSE before 29may2010
NULL # x$name <- x$name #[x$con$posZ]
else if ((is.na(x$name[1]))||(length(x$name)!=x$dim[2])) {
x$name <- sprintf(sprintf("%s%%0%dd"
,ifelse(is.na(x$name[1]),"Band ",x$name[1])
,nchar(length(1:x$con$bands))),1:x$con$bands)
# x$name <- character()
}
if (x$con$driver %in% c("ENVI","EGDAL"))
.write.hdr(x)
else if (x$con$driver=="RGDAL") {
.rgdal_prepare_con(x)
}
x$value <- NA
# class(x$value) <- ifelse(isCT,"ursaCategory","ursaNumeric")
class(x$value) <- ifelse(.is.colortable(x),"ursaCategory","ursaNumeric")
x
}
'.create.con' <- function(x,arglist)
{
if (.is.con(x$con))
{
if (is.null(arglist$datatype))
{
if ((TRUE)&&(!is.null(dim(x$value))))
{
nodata <- if (is.null(arglist$nodata)) x$con$nodata else arglist$nodata
if ((is.na(nodata))&&(anyNA(x$value)))
arglist$nodata <- nodata <- .optimal.nodata(x$value)
arglist$datatype <- .optimal.datatype(x$value,nodata)
if (!is.na(nodata)) {
if ((arglist$datatype==2)&&((nodata>32767)||(nodata<(-32768))))
arglist$nodata <- -32768L
else if ((arglist$datatype==12)&&((nodata>65535)||(nodata<0)))
arglist$nodata <- 65535L
else if ((arglist$datatype==1)&&((nodata>255)||(nodata<0)))
arglist$nodata <- 255L
else if ((arglist$datatype==11)&&((nodata>127)||(nodata<(-128))))
arglist$nodata <- -128L
}
}
else
arglist$datatype <- x$con$datatype ## release: only these line
}
if (is.null(arglist$nodata))
arglist$nodata <- x$con$nodata
if (is.null(arglist$connection))
arglist$connection <- "file"
}
else if (!is.null(dim(x$value)))
{
if (is.null(arglist$datatype))
arglist$datatype <- .optimal.datatype(x$value,x$con$nodata)
if (is.null(arglist$nodata)) {
arglist$nodata <- .make.nodata(arglist$datatype)
}
}
obj <- .raster.skeleton()
if (is.null(x))
{
obj$grid <- session_grid()
# obj$con <- envi=.con.skeleton()
obj$con <- .con.skeleton()
obj$con$driver <- arglist$implement ## added 20170124
# obj$con$driver <- driver ## removed 20170124
}
else
{
class(x) <- "ursaRaster (no generic)"
for (i in seq(along=x))
{
myname <- names(x)[i]
if (myname!="value")
obj[[myname]] <- x[[myname]]
}
}
con <- .make.con(obj,arglist)
con
}
'.make.con' <- function(obj,arglist)
{
if (!is.ursa(obj))
return(NULL)
if (missing(arglist))
arglist <- list()
grid <- obj$grid
if (!.is.con(obj$con))
con <- .con.skeleton()
else
con <- obj$con
ignore <- NULL #obj$nodata
fname <- NULL
datatype <- NULL
byteorder <- NULL
bands <- NULL
bandnames <- NULL
connectionName <- c("file","bzfile","gzfile","xzfile")
connection <- NULL
interleaveName <- c("bsq","bil","bip")
interleave <- NULL
implementName <- c("ENVI",ifelse(.isPackageInUse(),"EGDAL","RGDAL"))
implement <- NULL
driver <- NULL
proj <- NULL
compressed <- NULL
wkt <- NULL
ext <- NULL
for (i in seq(along=arglist))
{
var <- arglist[[i]]
if (is.null(var))
next
name <- names(arglist[i])
if ((is.null(fname))&&(is.character(var))&&
((!nchar(name))||(length(.grep("fname",name)))))
{
fname <- var
next
}
if ((is.null(connection))&&(is.character(var))&&
((!length(name))||(length(.grep("con",name)))))
{
if (length(.grep("^gz",var)))
connection <- "gzfile"
else if (length(.grep("^bz",var)))
connection <- "bzfile"
else if (length(.grep("^xz",var)))
connection <- "xzfile"
else if (length(.grep("^file",var)))
connection <- "file"
if (!is.null(connection))
next
}
if ((is.null(interleave))&&(is.character(var))&&
((!length(name))||(length(.grep("interleave",name)))))
{
ind <- .grep(var,interleaveName)
if (length(ind))
{
interleave <- interleaveName[ind]
next
}
}
if ((is.null(implement))&&(is.character(var))&&
((!length(name))||(length(.grep("implement",name)))))
{
ind <- .grep(var,implementName)
if (length(ind))
{
implement <- implementName[ind]
next
}
}
if ((is.null(datatype))&&
((!length(name))||(length(.grep("datatype",name)))))
{
if (is.numeric(var))
{
if (as.integer(var) %in% c(1,2,3,4,5,12,13,11))
{
datatype <- var
next
}
}
else if (is.character(var))
{
datatype <- switch(var,byte=1L,integer=2L,real=4L,float=4L
,Byte=1L,UInt8=1L,Int8=11
,Int16=2L,UInt16=12,UInt32=13,Int32=3
,Float32=4L,Float64=5L
,NULL)
next
}
}
if ((is.null(ignore))&&(!is.ursa(var))&&((is.numeric(var))||(anyNA(var)))&&
(length(name))&&(length(.grep("(^bg$|nodata|ignore)",name))))
{
ignore <- var
next
}
if ((is.null(byteorder))&&(is.numeric(var))&&
((length(name))&&(length(.grep("byteorder",name)))))
{
byteorder <- as.integer(var)
if (byteorder)
byteorder <- 1L
next
}
# if ((is.null(bands))&&(is.numeric(var))&&((length(name))
# &&((length(.grep(name,"bands")))||(length(.grep(name,"nband")))||
# (length(.grep(name,"layers")))||(length(.grep(name,"nlayer"))))))
if ((is.null(bands))&&(is.numeric(var))&&(length(name))&&
(length(.grep("(layers|bands|nlayer|nband|length)",name))))
{
bands <- as.integer(var)
next
}
# if ((is.null(bandnames))&&(is.character(var))&&((length(name))
# &&((length(.grep(name,"bandname")))||
# (length(.grep(name,"layername")))||
# (length(.grep(name,"name"))))))
if ((is.null(bandnames))&&(is.character(var))&&(length(name))&&
(length(.grep("(bandname|layername|names)",name))))
{
bandnames <- var
next
}
if ((is.null(driver))&&(is.character(var))&&(length(name))&&
(length(.grep("driver",name))))
{
driver <- var
next
}
if ((is.null(compressed))&&((is.integer(var) || is.logical(var)))&&
((length(name))&&(length(.grep("compress",name)))))
{
compressed <- as.logical(var)
next
}
if ((is.null(wkt))&&((is.integer(var) || is.logical(var)))&&
((length(name))&&(length(.grep("wkt",name)))))
{
wkt <- as.logical(var)
next
}
if ((is.null(ext))&&(is.character(var))&&
((length(name))&&(length(.grep("ext",name)))))
{
ext <- var
next
}
}
if ((is.null(bands))&&(length(bandnames)))
bands <- length(bandnames)
if (is.null(fname))
{
fname <- .maketmp()
message(paste("Filename is assigned automatically:",fname))
}
myDir <- .dirname(fname)
if (!file.exists(myDir)) {
dir.create(myDir,recursive=TRUE)
}
con$fname <- file.path(chartr("\\","/",normalizePath(.dirname(fname)))
,.basename(fname))
if (!is.null(connection))
con$connection <- connection
else
connection <- "file"
if (is.na(con$connection))
con$connection <- connection
if (!is.null(datatype))
con$datatype <- datatype
else
datatype <- 4L
if (is.na(con$datatype))
con$datatype <- datatype
if (!is.null(ignore))
con$nodata <- ignore
else
ignore <- .make.nodata(datatype)
if (is.na(con$nodata))
con$nodata <- ignore
if (!is.null(interleave))
con$interleave <- interleave
else
interleave <- interleaveName[1]
if (is.na(con$interleave))
con$interleave <- interleave
if (!is.null(implement))
con$driver <- implement
else
implement <- implementName[1]
if (is.na(con$driver))
con$driver <- implement
if (!is.null(byteorder))
con$byteorder <- byteorder
else
byteorder <- 0L
if (is.na(con$byteorder))
con$byteorder <- byteorder
if ((con$driver %in% c("ENVI","EGDAL"))&&(TRUE)) ## forced to compress
{
if (!(is.null(compressed))) {
if (!is.na(compressed))
con$compress <- as.integer(compressed)
else
con$compress <- 1L
}
else
con$compress <- 1L
}
else
{
if (!is.null(compressed))
con$compress <- as.integer(compressed)
else
compressed <- 0L
if (is.na(con$compress))
con$compress <- compressed
con$compress <- abs(con$compress)
}
con$compress <- abs(con$compress)
if (!is.null(wkt))
con$wkt <- as.logical(wkt)
else
wkt <- FALSE
if (is.na(con$wkt))
con$wkt <- wkt
if (is.null(ext))
ext <- switch(con$driver,ENVI=".envi","")
else
{
ext <- if (nchar(ext)) paste0(".",ext) else ""
con$compress=0L
}
if (is.na(con$wkt))
con$wkt <- wkt
if (!is.null(bands))
con$bands <- bands
else if (length(obj$value))
bands <- obj$dim[2]
else
bands <- 1L
if (is.na(con$bands))
con$bands <- bands
con$mode <- if (con$datatype %in% c(1L,2L,3L,11L,12L,13L)) "integer"
else "numeric"
con$sizeof <- if (con$datatype %in% c(1L,11L)) 1L
else if (con$datatype %in% c(2L,12L)) 2L
else if (con$datatype %in% c(3L,13L,4L)) 4L else 8L
con$signed <- if (con$datatype %in% c(1L,12L,13L)) FALSE else TRUE
con$endian <- "little"
con$swap <- as.integer((.Platform$endian=="big")&(!con$byteorder)|
(.Platform$endian=="little")&(con$byteorder))
# endian != .Platform$endian
if (FALSE)
NULL
else if (length(.grep("\\.gz$",con$fname)))
con$connection <- "gzfile"
else if (length(.grep("\\.bz2$",con$fname)))
con$connection <- "bzfile"
else if (length(.grep("\\.xz$",con$fname)))
con$connection <- "xzfile"
else if (length(.grep("\\.(envi|bin|img|dat|tif|tiff|hfa)$",con$fname)))
{
con$connection <- "file"
con$compress <- 0L
}
##~ else if (length(.grep("\\.bin$",con$fname)))
##~ {
##~ con$connection <- "file"
##~ con$compress <- 0L
##~ }
##~ else if (length(.grep("\\.envi$",con$fname)))
##~ {
##~ con$connection <- "file"
##~ con$compress <- 0L
##~ }
##~ else if (length(.grep("\\.img$",con$fname)))
##~ {
##~ con$connection <- "file"
##~ con$compress <- 0L
##~ }
else if (length(.grep("\\.$",con$fname)))
{
con$fname <- substr(con$fname,1,nchar(con$fname)-1L)
con$connection <- "file"
con$compress <- 0L
}
else
con$fname <- sprintf("%s%s",con$fname,switch(con$connection
,file=ext,gzfile=".gz",bzfile=".bz2",xzfile=".xz",""))
if (!(con$connection %in% c("file")))
con$compress <- 0L
con$lines <- grid$rows
con$samples <- grid$columns
# con$offset <- as.integer(5*4+8*8+sum(nchar(obj$name))+length(obj$name)+
# nchar(obj$grid$crs)+1+
# 4*obj$dim[2]+con$sizeof)
if ((length(con$offset)>1)||(is.na(con$offset)))
con$offset <- 0L
if (is.na(con$bands))
{
if (!is.null(bandnames))
con$bands <- length(bandnames)
else
con$bands <- 1L
}
if ((!is.na(obj$dim[2]))&&(is.null(bands)))
con$bands <- obj$dim[2]
con$interleave <- with(con,switch(interleave,spatial="bsq",temporal="bil"
,interleave))
if (con$driver %in% c("ENVI","EGDAL")) {
if (con$driver=="ENVI")
con$handle <- do.call(con$connection,list(con$fname,"w+b"))
else {
ftemp <- .maketmp(1)
con$fname <- c(con$fname,ftemp)
con$handle <- do.call(con$connection,list(ftemp,"w+b"))
}
cl <- class(con$handle)
if (("bzfile" %in% cl)||("xzfile" %in% cl)||("gzfile" %in% cl)) ## gz?
con$seek <- FALSE
else
con$seek <- TRUE
rm(cl)
}
else if (con$driver %in% c("EGDAL","RGDAL")) {
if (.lgrep("\\.tif(f)*$",fname))
driver <- "GTiff"
else if (.lgrep("\\.img$",fname))
driver <- "HFA" # https://gdal.org/frmt_hfa.html
else if (.lgrep("\\.png$",fname))
driver <- "PNG"
else if (.lgrep("\\.jp(e)*g$",fname))
driver <- "JPEG"
else if (.lgrep("\\.bmp$",fname))
driver <- "BMP"
# else if (.lgrep("\\.sdat$",fname))
# driver <- "SAGA"
if (is.null(driver))
driver <- "ENVI"
dtName <- switch(as.character(datatype)
,'1'="Byte",'2'="Int16",'4'="Float32"
,'11'="Int8",'12'="UInt16",'13'="UInt32",'3'="Int32"
,'5'="Float64",stop("cannot recognize datatype"))
nb <- if (is.na(con$posZ[1])) con$bands else length(con$posZ)
if (con$driver %in% "RGDAL") {
try(con$handle <- methods::new("GDALTransientDataset"
,methods::new("GDALDriver",driver)
,con$lines,con$samples,nb,dtName))
if (!inherits(con$handle,"GDALTransientDataset"))
return(NA)
con$seek <- FALSE
}
else { ## via ENVI -> GDAL after closing
ftemp <- .maketmp(1)
con$fname <- c(con$fname,ftemp)
con$handle <- do.call(con$connection,list(ftemp,"w+b"))
con$seek <- TRUE
}
}
if ("file" %in% class(con$handle))
{
if (file.exists(oldpack <- paste0(con$fname[1],".gz")))
file.remove(oldpack)
if (file.exists(oldpack <- paste0(con$fname[1],".bz2")))
file.remove(oldpack)
if (file.exists(oldpack <- paste0(con$fname[1],".xz")))
file.remove(oldpack)
}
if (file.exists(aux <- paste0(con$fname[1],".aux.xml")))
file.remove(aux)
if (con$driver %in% c("ENVI","EGDAL")) {
if (((con$interleave %in% c("bil","bsq"))&&(con$seek))||
(con$connection %in% c("gzfile"))) #&&(!is.matrix(obj$value)))
{
nb <- ifelse(is.na(con$posZ[1]),con$bands,length(con$posZ))
seek(con,origin="start"
,where=with(con,(lines*samples*nb-1)*sizeof+offset),rw="w")
val <- 0L
storage.mode(val) <- con$mode
with(con,writeBin(val,size=sizeof,endian=endian,handle))
}
}
con$indexC <- NA
con$indexR <- NA
if ((!is.na(con$posZ[1]))&&(!is.na(con$indexZ[1])))
con$posZ <- con$indexZ[con$posZ]
con$indexZ <- NA
storage.mode(con$nodata) <- con$mode
storage.mode(con$datatype) <- "integer"
# con$driver <- "ENVI"
con$swap <- as.integer((.Platform$endian=="big")&(!con$byteorder)|
(.Platform$endian=="little")&(con$byteorder))
con
}
'.optimal.datatype' <- function(x,nodata=NULL)
{
if (is.ursa(x))
x <- x$value
if (TRUE) {
isInt <- is.integer(x)
n <- as.integer(prod(dim(x)))+as.integer(!is.null(nodata))
if (isInt) {
res <- .Cursa(C_optimalDatatypeInt,x=c(as.integer(x),as.integer(nodata))
,n=n,res=integer(1),NAOK=TRUE)$res
}
else {
res <- .Cursa(C_optimalDatatypeDouble,x=c(as.numeric(x),as.numeric(nodata))
,n=n,res=integer(1),NAOK=TRUE)$res
}
# str(res)
# print(summary(x));q()
return(res)
}
# isInteger <- is.integer(x)
dimx <- dim(x)
lend <- length(dimx)
if (!(length(dimx) %in% c(2,3)))
return(NULL)
# if (!isInteger)
# datatype <- 4
bands <- if (lend==3) dimx[3] else if (lend==2) dimx[1]
if (is.na(bands))
return(NULL)
datatypeout <- rep(NA_integer_,bands)
tmp <- if (lend==3) array(NA,dim=dimx[-lend]) else array(NA,dim=dimx[2])
for (i in seq(bands))
{
tmp[] <- if (lend==3) x[,,i] else if (lend==2) x[i,]
tmp[is.na(tmp)] <- 0
if ((all(tmp %in% 0:255))||(all(tmp %in% -128:127)))
datatypeout[i] <- 1L
else if (all(tmp %in% -32768:32767))
datatypeout[i] <- 2L
else if (all(tmp %in% 0:65536))
datatypeout[i] <- 12L
else
datatypeout[i] <- 4L
}
rm(tmp)
datatype <- unique(datatypeout)
if (length(datatype)>1)
{
if (4L %in% datatype)
datatype <- 4L
else if ((1L %in% datatype)&&(2L %in% datatype)&&(!(12L %in% datatype)))
datatype <- 2L
else if ((1L %in% datatype)&&(12L %in% datatype)&&(!(2L %in% datatype)))
datatype <- 12L
else
datatype <- 4L
}
datatype
}
'.write.hdr' <- function(x,clear=TRUE)
{
if (!is.ursa(x))
return(NULL)
grid <- x$grid
con <- x$con
if (!.is.grid(grid))
return(NULL)
if (!.is.con(con))
return(NULL)
if (is.na(con$interleave))
con$interleave <- "bsq"
if (is.na(con$samples))
con$samples <- grid$columns
if (is.na(con$lines))
con$lines <- grid$rows
if (is.na(con$bands))
con$bands <- 1L
if (is.na(con$byteorder))
con$byteorder <- 0L
if (is.na(con$datatype))
con$datatype <- 1L
indF <- ifelse(length(con$fname)>1,2L,1L)
fname <- .gsub("\\.(bin|bingz|envi|envigz|img|dat|gz|bz2|xz|unpacked(.*)~)$",""
,con$fname[indF])
if (is.na(fname))
{
fname <- x$name[1]
con$fname <- fname
}
myname <- paste(fname,".hdr",sep="")
if (clear)
{
list1 <- .dir(path=dirname(fname)
,pattern=sprintf("^%s($|\\.(envi|envigz|bin|bingz|dat|img|gz|bz2|xz|.*aux\\.xml)$)"
,basename(fname)),full.names=TRUE)
if (length(ind <- .grep(sprintf("^%s$",basename(con$fname[indF])),basename(list1))))
list1 <- list1[-ind]
file.remove(list1[which(!file.info(list1)$isdir)])
}
Fout <- file(myname,"wt")
writeLines(sprintf("%s","ENVI"),Fout)
writeLines(sprintf("description = {%s}",con$fname[indF]),Fout)
writeLines(sprintf("samples = %d",con$samples),Fout)
writeLines(sprintf("lines = %d",con$lines),Fout)
if (is.na(con$posZ[1]))
writeLines(sprintf("bands = %d",con$bands),Fout)
else
writeLines(sprintf("bands = %d",length(con$posZ)),Fout)
writeLines(sprintf("header offset = %d",0),Fout)
writeLines(sprintf("file type = %s","ENVI Standard"),Fout)
writeLines(sprintf("data type = %d",con$datatype),Fout)
writeLines(sprintf("interleave = %s",con$interleave),Fout)
if (FALSE)
writeLines(sprintf("sensor type = %s","Unknown"),Fout)
writeLines(sprintf("byte order = %d",con$byteorder),Fout)
if ((length(unique(con$nodata))==1)&&(!is.na(con$nodata)))
writeLines(sprintf("data ignore value = %s",con$nodata[1]),Fout)
projection_name <- "ursa package"
projection_ellipse <- "unknown"
projection_units <- "units=meters"
projection_info <- ""
proj4 <- grid$crs[1]
p <- NULL
if (nchar(proj4))
{
pr <- unlist(strsplit(proj4,"\\s+"))
projection_ellipse <- "unknown"
projection_units <- "units=meters"
projection_name <- "R"
semi <- NULL
if (FALSE)
semi <- c("6378137","6356752.3") ## low-precision
if (.lgrep("\\+datum=WGS84",pr))
semi <- c("6378137","6356752.314245179")
else if (.lgrep("\\+ellps=WGS84",pr))
semi <- c("6378137","6356752.314245179")
else if (.lgrep("\\+ellps=krass",pr))
semi <- c("6378245","6356863.018773047")
else
semi <- c(.gsub2("\\s\\+a=(\\S+)\\s","\\1",proj4)
,.gsub2("\\s\\+b=(\\S+)\\s","\\1",proj4))
if (.lgrep("\\+proj=laea",pr))
{
p <- c(p,"11",semi)
p <- c(p,.gsub2("\\s\\+lat_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lon_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+x_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+y_0=(\\S+)\\s","\\1",proj4))
projection_name <- "Lambert Azimuthal Equal Area"
}
else if (.lgrep("\\+proj=stere",pr))
{
if (.lgrep("\\+lat_ts",proj4)) {
p <- c(p,"31",semi)
p <- c(p,.gsub2("\\s\\+lat_ts=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lon_0=(\\S+)(\\s|$)","\\1",proj4))
if(.lgrep("\\s\\+x_0=(\\S+)(\\s|$)",proj4))
p <- c(p,.gsub2("\\s\\+x_0=(\\S+)(\\s|$)","\\1",proj4))
else
p <- c(p,"0")
if(.lgrep("\\s\\+y_0=(\\S+)(\\s|$)",proj4))
p <- c(p,.gsub2("\\s\\+y_0=(\\S+)(\\s|$)","\\1",proj4))
else
p <- c(p,"0")
projection_name <- "Polar Stereographic"
}
else {
p <- c(p,"7",semi)
p <- c(p,.gsub2("\\s\\+lat_0=(\\S+)(\\s|$)","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lon_0=(\\S+)(\\s|$)","\\1",proj4))
if(.lgrep("\\s\\+x_0=(\\S+)(\\s|$)",proj4))
p <- c(p,.gsub2("\\s\\+x_0=(\\S+)(\\s|$)","\\1",proj4))
else
p <- c(p,"0")
if(.lgrep("\\s\\+y_0=(\\S+)(\\s|$)",proj4))
p <- c(p,.gsub2("\\s\\+y_0=(\\S+)(\\s|$)","\\1",proj4))
else
p <- c(p,"0")
projection_name <- "Stereographic (ellipsoid)"
}
}
else if (.lgrep("\\+proj=tmerc",pr))
{
p <- c(p,"3",semi)
p <- c(p,.gsub2("\\s\\+lat_ts=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lon_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+x_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+y_0=(\\S+)\\s","\\1",proj4))
projection_name <- "Transverse Mercator"
}
else if (.lgrep("\\+proj=lcc",pr))
{
p <- c(p,"4",semi)
p <- c(p,.gsub2("\\s\\+lat_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lon_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+x_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+y_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lat_1=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lat_2=(\\S+)\\s","\\1",proj4))
projection_name <- "Lambert Conformal Conic"
}
else if (.lgrep("\\+proj=aea",pr))
{
p <- c(p,"9",semi)
p <- c(p,.gsub2("\\s\\+lat_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lon_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+x_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+y_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lat_1=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lat_2=(\\S+)\\s","\\1",proj4))
projection_name <- "Albers Conical Equal Area"
}
else if (.lgrep("\\+proj=zzzcea",pr))
{
p <- c(p,"99",semi)
p <- c(p,.gsub2("\\s\\+lon_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+lat_ts=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+x_0=(\\S+)\\s","\\1",proj4))
p <- c(p,.gsub2("\\s\\+y_0=(\\S+)\\s","\\1",proj4))
projection_name <- "Equal Area Cylindrical"
}
if ((FALSE)||(!nchar(projection_name)))
{
print(pr)
print(projection_info)
stop("NEW PROJECTION")
}
if (!is.null(p))
{
p <- c(p,projection_name,projection_units)
projection_info <- paste0("{",paste(p,collapse=", "),"}")
}
}
tempx <- with(grid,minx+0*0.5*resx)
tempy <- with(grid,grid$maxy-0*0.5*resy)
map_info <- c(projection_name,1,1,tempx,tempy,grid$resx,grid$resy
,projection_ellipse,projection_units)
map_info <- paste0("{",paste(map_info,collapse=", "),"}")
writeLines(sprintf("map info = %s",map_info),Fout)
if ((TRUE)&&(nchar(Sys.which("gdalsrsinfo")))) ## forced 'gdalsrsinfo'/showWKT
projection_info <- ""
if (projection_info!="")
writeLines(sprintf("projection info = %s",projection_info),Fout)
wkt <- NULL
if ((TRUE)&&(!is.null(p))) ## help QGIS to recognize EPSG
{
.epsg <- NA
if ((p[1]=="31")&&
(as.numeric(p[2])==6378137)&&(as.numeric(p[3])==6356752.314245179)&&
(as.numeric(p[4])==70)&&(as.numeric(p[5])==-45)) {
.epsg <- "3408"
wkt <- paste0("PROJCS[\"WGS_84_NSIDC_Sea_Ice_Polar_Stereographic_North\""
,",GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\""
,",SPHEROID[\"WGS_1984\",6378137,298.257223563]]"
,",PRIMEM[\"Greenwich\",0]"
,",UNIT[\"Degree\",0.017453292519943295]]"
,",PROJECTION[\"Stereographic_North_Pole\"]"
,",PARAMETER[\"standard_parallel_1\",70]"
,",PARAMETER[\"central_meridian\",-45]"
,",PARAMETER[\"false_easting\",0]"
,",PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]}")
}
else if ((p[1]=="36")&&
(as.numeric(p[2])==6378137)&&
(as.numeric(p[3])==6356752.314245179)&&
(as.numeric(p[4])==90))
{
if (as.numeric(p[5])==180)
{
.epsg <- "3571"
wkt <- paste0("PROJCS[\"WGS_84_North_Pole_LAEA_Bering_Sea\""
,",GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\""
,",SPHEROID[\"WGS_1984\",6378137,298.257223563]]"
,",PRIMEM[\"Greenwich\",0]"
,",UNIT[\"Degree\",0.017453292519943295]]"
,",PROJECTION[\"Lambert_Azimuthal_Equal_Area\"]"
,",PARAMETER[\"latitude_of_origin\",90]"
,",PARAMETER[\"central_meridian\",180]"
,",PARAMETER[\"false_easting\",0]"
,",PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]")
}
else if (as.numeric(p[5])==-150) {
.epsg <- "3572"
wkt <- paste0("PROJCS[\"WGS_84_North_Pole_LAEA_Alaska\""
,",GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\""
,",SPHEROID[\"WGS_1984\",6378137,298.257223563]]"
,",PRIMEM[\"Greenwich\",0]"
,",UNIT[\"Degree\",0.017453292519943295]]"
,",PROJECTION[\"Lambert_Azimuthal_Equal_Area\"]"
,",PARAMETER[\"latitude_of_origin\",90]"
,",PARAMETER[\"central_meridian\",-150]"
,",PARAMETER[\"false_easting\",0]"
,",PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]")
}
else if (as.numeric(p[5])==90) {
.epsg <- "3576"
wkt <- paste0("PROJCS[\"WGS_84_North_Pole_LAEA_Russia\""
,",GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\""
,",SPHEROID[\"WGS_1984\",6378137,298.257223563]]"
,",PRIMEM[\"Greenwich\",0]"
,",UNIT[\"Degree\",0.017453292519943295]]"
,",PROJECTION[\"Lambert_Azimuthal_Equal_Area\"]"
,",PARAMETER[\"latitude_of_origin\",90]"
,",PARAMETER[\"central_meridian\",90]"
,",PARAMETER[\"false_easting\",0]"
,",PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]")
}
}
# print(.epsg)
rm(.epsg)
}
if ((is.null(wkt))&&(projection_info=="")&&(nchar(proj4)))
{
lverbose <- FALSE
if (lverbose)
.elapsedTime("proj4 -> wkt start")
# loadNamespace("sf") ## for development
# (!("package:rgdal" %in% search()))) {
if ((nchar(Sys.which("gdalsrsinfo")))&&
(!(any(c("rgdal","sf") %in% loadedNamespaces())))) {
if (lverbose)
message("'gdalsrsinfo' engine (write)")
if (FALSE)
wkt <- system2("gdalsrsinfo",c("-o wkt",.dQuote(proj4))
,stdout=TRUE)
else {
wktout <- .maketmp(ext=".wkt~")
# shell(paste("gdalsrsinfo -o wkt",paste0("\"",proj4,"\""),"1>",wktout))
# 20170319 dQuote() returns non-symmetrical quotes in interactive()
system2("gdalsrsinfo",c("-o wkt_esri",.dQuote(proj4))
,stdout=wktout,stderr=FALSE)
wkt <- readLines(wktout,warn=FALSE)
file.remove(wktout)
}
}
else {
if (!("sf" %in% loadedNamespaces())) {
if (.isPackageInUse()) {
if (lverbose)
message("forced 'sf' loading")
requireNamespace("sf") ## rgdal retired
}
}
if (!("sf" %in% loadedNamespaces())) {
if (lverbose)
message("'rgdal' engine")
if (!.try(wkt <- .rgdal_showWKT(proj4,morphToESRI=TRUE)))
wkt <- NULL
}
else { ## 'sf' in namespace; 'OGC_WKT' ONLY.
if (lverbose)
message("'sf' engine")
if (!.try(wkt <- {
if (utils::packageVersion("sf")<"0.9")
ret <- sf::st_as_text(sf::st_crs(proj4),EWKT=TRUE)
else
ret <- sf::st_crs(proj4)$Wkt
ret
}))
# if (!.try(wkt <- sf::st_as_text(sf::st_crs(proj4),EWKT=TRUE)))
# wkt <- NULL
# print(proj4)
# message(wkt)
if (!TRUE) { ## 20191216 patch for EXTENSION["PROJ4","+proj=......."]
wkt <- gsub(",EXTENSION\\[\"PROJ4\".+\\]","]",wkt)
}
}
}
if (lverbose)
.elapsedTime("proj4 -> wkt finish")
}
if (!is.null(wkt)) {
if (length(wkt))
wkt <- paste(gsub("(^\\s+|\\s+$)","",wkt),collapse="")
writeLines(sprintf("coordinate system string = {%s}",wkt),Fout)
}
if ((is.character(x$name))&&(sum(nchar(x$name))>0))
{
if (.lgrep(",",x$name)) {
metaname <- paste0(.gsub("\\.(gz|bz2|xz)$","",x$con$fname[indF]),".aux.xml")
Fmeta <- file(metaname,"wt")
writeLines("<PAMDataset>",Fmeta)
writeLines(" <Metadata>",Fmeta)
for (i in seq_along(x$name))
writeLines(paste0(" <MDI key=",.dQuote(paste0("Band_",i))
,">",x$name[i],"</MDI>"),Fmeta)
writeLines(" </Metadata>",Fmeta)
writeLines("</PAMDataset>",Fmeta)
close(Fmeta)
}
else {
if (all(!is.na(x$con$posZ)))
bandnames <- paste(x$name[x$con$posZ],sep="",collapse=",\n\t")
else
bandnames <- paste(x$name,sep="",collapse=",\n\t")
writeLines(sprintf("band names = {\n\t%s}",bandnames),Fout)
}
}
if (length(x$colortable))
{
ct <- x$colortable
writeLines(sprintf("classes = %d",length(ct)),Fout)
isCol <- TRUE
if (is.numeric(ct))
val <- as.numeric(ct)
else {
val <- as.character(ct)
if (all(is.na(val)))
isCol <- FALSE
else
val <- c(col2rgb(as.character(val)))
}
if (isCol) {
val <- paste(val,sep="",collapse=",")
writeLines(sprintf("class lookup = {%s}",val),Fout)
}
if (!is.null(names(ct))) {
if (TRUE) ## if '\t' not recoginizing
{
val <- .gsub(",",";",names(x$colortable))
val <- paste(val,sep="",collapse=",")
writeLines(sprintf("class names = {%s}",val),Fout)
}
else
{
val <- paste(names(x$colortable),sep="",collapse=",\n\t")
writeLines(sprintf("class names = {\n\t%s}",val),Fout)
}
}
}
close(Fout)
invisible(chartr("\\","/",normalizePath(myname)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.