R/offline.R

In rgplates: R Interface for the GPlates Web Service and Desktop Application

###########################################################################
# Offline interface to Gplates


#Mac examples
#	x<- fetch("pared", "public")[,c("longitude", "latitude")]
#	mo <- fetch("paleomap", "model")
#	reconstruct(x, age=10, model=mo, verbose=TRUE)
#reconstruct(x, age=10, model=mo, verbose=TRUE, path.gplates="/Users/Nick/Downloads/GPlates-2.2.0/gplates.app/Contents/MacOS/gplates")

reconstructGPlates <- function(x, age, model, path.gplates=NULL,dir=NULL, verbose=FALSE, cleanup=TRUE, plateperiod=FALSE, gmeta=FALSE, partitioning="static_polygons", check=TRUE){
	# shortcut - will be put in a matrix automatically
	if(is.na(age)) return(NA)
	if(!inherits(model, "platemodel")) stop("You need a GPlates tectonic object for this method.")
	# 0. Process the attributes
	if(inherits(model@features, "data.frame")){
		if(check){
			# check for the ages!
			if(is.character(x)){
				checkThis <- x
			#fall through, use the partitioning!
			}else{
				checkThis <- partitioning
			}
			# the part corresponding to this
			if(!any(checkThis==rownames(model@features))) stop(paste0("The feature collection '",checkThis,"' is not found."))
			thisPart <-model@features[ checkThis, ] 

			# get the ages
			young <- thisPart[ , "to"]
			old <- thisPart[ , "from"]
			if(!(age>=young && age <=old) ) stop(paste0("Invalid age. '", checkThis, "' is valid between ", old, " Ma and ", young, " Ma."))
		}

		# omit the age info, use characters as features
		nam <- rownames(model@features)
		model@features <- model@features[, "feature_collection"]
		names(model@features) <- nam
	}	
	
    # 1. FIND GPlates
		# A. get operating system
		os <- getOS()

		# assume a non-retarded operating system...
		win <- FALSE

		# B. should the program look for a default path for gplates?
		if(is.null(path.gplates)){
			# depending on the os
			if(os=="linux"){
				# the GPlates executable itself
				gplatesExecutable <- "gplates"

				# what the user would have entered
				path.gplates <- gplatesExecutable

				# leave the model intact in the namespace (easier debug)
				rotation <- model@rotation
				plateFeatures<- model@features

				# separator character between directories
				dirSep <- "/"

			}
			if(os=="windows"){
				win <- TRUE
				# 1. find GPLATES exectutable if possible
				# directory and executable
				gplatesPaths <- winDefaultGPlates()
				#path to executable
				path.gplates <- paste(gplatesPaths, collapse="/")
				# system call to executable
				gplatesExecutable <- paste("\"", path.gplates, "\"", sep="")

				# 2. replace model paths with \\
				rotation <- gsub("/","\\\\", model@rotation)
				plateFeatures<- gsub("/","\\\\", model@features)

				# characters to include directory 
				dirSep <- "\\\\"

			}
			if(os=="osx"){
				# the GPlates executable itself
				# default
				# gplatesExecutable <- "/Applications/GPlates-2.2.0/gplates.app/Contents/MacOS/gplates"
				gplatesPaths <- macDefaultGPlates()
				gplatesExecutable <- paste(gplatesPaths, collapse="/")
				
				# what the user would have entered
				path.gplates <- gplatesExecutable

				# leave the model intact in the namespace (easier debug)
				rotation <- model@rotation
				plateFeatures <- model@features

				# separator character between directories
				dirSep <- "/"
				
			}

		# look for given path
		}else{
			# separate to form a length 2 vector
			gplatesExecutable <- path.gplates
			
			# leave the model intact in the namespace (easier debug)
			rotation <- model@rotation
			plateFeatures <- model@features

			# separator character between directories
			dirSep <- "/"
			
			# windows needs special treatment
			if(os=="windows"){
				win <- TRUE
			
				# system call to executable
				gplatesExecutable <- paste("\"", path.gplates, "\"", sep="")

				# 2. replace model paths with \\
				rotation <- gsub("/","\\\\", model@rotation)
				plateFeatures <- gsub("/","\\\\", model@features)

				# characters to include directory 
				dirSep <- "\\\\"

			}
			
		}

		# C. one moretest whether gplates was detected or not
			gpTest <- testGPlates(gplatesExecutable, verbose=verbose)

			# if gplates is not present:
			if(!gpTest) stop(paste("The GPlates executable\n	\"", path.gplates,"\"\nwas not found.", sep=""))

	# 2. Setup reconstruction environment
		# folder where files will be executed
		if(is.null(dir)) tempd <- tempdir() else tempd <- dir

		# copy all model feature to working directory
		sources <- plateFeatures
		plateFeatures <-file.path(tempd, fileFromPath(sources, win=win)) 
		# bug fix for windows
		if (win) plateFeatures <- gsub("/","\\\\", plateFeatures)

		names(plateFeatures) <- names(sources)
		results <- file.copy(sources, plateFeatures)

		# prepare x
		# create a SpatialPointsDataFrame from long-lat matrix

		originals <- NULL 
		if(inherits(x, "matrix")){
			originals <- x
			x <- as.data.frame(x)
		}
			
		# data.frame
		if(inherits(x, "data.frame") & !(inherits(x, "sf") | inherits(x, "Spatial"))){
			# the original
			if(is.null(originals)) originals <- x
			# enforce column names
			colnames(x) <- c("lng", "lat")
			x$ID <- paste0("a", 1:nrow(x))
			theID <- x$ID

			# transform to sf
			x <- sf::st_as_sf(x, coords=c("lng", "lat"))
		}
		
	
		# spatial original, enforce sf!
		if(inherits(x, "Spatial")){
			# conserve the originals to know what to have
			originals <- x
			x <- sf::st_as_sf(x)
		}
	
		# for non -spatial/sf this is null
		originalCRS <- NULL
		# enforce projection!
		# if originally an sf!
		if(inherits(x, "sf")){
			if(is.null(originals)) originals <- x
			# the original CRS
			originalCRS <- sf::st_crs(x)
			if(!is.na(originalCRS)){
				xTransform <- sf::st_transform(x, "EPSG:4326")
			}else{
				sf::st_crs(x) <- "EPSG:4326"
				xTransform <- x
			}
		}


		# in case stat
		if(!is.character(x)){
		# write 'x' as a shapefile
			layer<- paste(randomString(length=3), age, sep="_")
			if(verbose) message(paste("Exported data identified as ", layer))
			pathToFileNoEXT <- paste(tempd, "/", layer,sep="")
			if (win) pathToFileNoEXT <- gsub("/","\\\\", pathToFileNoEXT)
			if(verbose) message("Exporting 'x' as a shapefile.")
#			rgdal::writeOGR(xTransform, dsn=paste(pathToFileNoEXT, ".shp", sep=""), layer=layer, driver="ESRI Shapefile")
			sf::st_write(xTransform, dsn=paste(pathToFileNoEXT, ".shp", sep=""),
				layer=layer, driver="ESRI Shapefile", quiet=!verbose)
			# select partitioning feature
			platePolygons <- plateFeatures[partitioning]


		}else{
		# feature to reconstruct is the static polygons
			if(length(x)!=1) stop("Only one feature collection can be reconstructed with this method.")
			# look for x in the feature set
			validFeature <- any(x==names(plateFeatures))
			if(validFeature){
				# use original one - even for windows.
				pathToFileNoEXT <- gsub(".gpml", "",plateFeatures[x])
				if (win) pathToFileNoEXT <- gsub("/","\\\\", pathToFileNoEXT)
			}else{
				stop("The requested feature collection is not part of the model. ")
			}
		}

		# inheritance of appearance and disappearance dates
		if(plateperiod){
			pPer <- 1
		}else{
			pPer <- 0
		}
		
	# 3. Execute GPlates commands
		# convert to gpml
		if(!is.character(x)){
			if(verbose) message("Converting shapefile to .gpml.")
			conversion <- paste(gplatesExecutable, " convert-file-format -l \"",pathToFileNoEXT,".shp\" -e gpml", sep="")
			system(conversion, ignore.stdout=!verbose,ignore.stderr=!verbose)
		}
		# do the plate assignment
		if(!is.character(x)){
			if(verbose) message("Assigning plate IDs to .gpml file.")
			assignment <- paste(gplatesExecutable, " assign-plate-ids -e ",pPer," -p \"", platePolygons, "\" -l \"",pathToFileNoEXT,".gpml\"", sep="")
			system(assignment, ignore.stdout=!verbose,ignore.stderr=!verbose)
		}
		# do reconstruction
		if(!is.character(x)) if(verbose) message("Reconstructing coordinates.")
		if(is.character(x)) if(validFeature) if(verbose) message(paste0("Reconstructing '",  x, "'."))
			reconstruction <- paste(gplatesExecutable, " reconstruct -l \"",pathToFileNoEXT,".gpml\" -r \"", 
					rotation, "\" -e shapefile -t ", age, " -o \"", pathToFileNoEXT,"_reconstructed\" -w 1", sep="") 
			system(reconstruction, ignore.stdout=!verbose,ignore.stderr=!verbose)

	# 4. Processing output
		# Was the output multiple?
		multiple <- FALSE

		# reading coordinates
		if(verbose) message("Reading reconstructed geometries.")

		# the single file
		targetSingle <- paste(pathToFileNoEXT,"_reconstructed.shx",	sep="")
		targetSingleNoPath <- fileFromPath(targetSingle, win=win)

		# produced directory? 
		targetDir<- paste(pathToFileNoEXT,"_reconstructed",	sep="")
		targetDirNoPath <- fileFromPath(targetDir, win=win)

		# is this a single file? 
		allFiles <- list.files(tempd)

		# assume that nothing was calculated. The formatting will be taken care of by the front-end functions
		rotated <- NA
	
		# is the target single file created?
		if(any(allFiles==targetSingleNoPath)){
			if(verbose) message("Found single output geometry files.")
			rotated <- sf::st_read(targetSingle, quiet=!verbose)
		}

		# did GPlates produce a whole directory? 
		if(any(allFiles==targetDirNoPath)){
			# Was an output multiple?
			multiple <- TRUE

			if(verbose) message("Found multiple output geometry files.")
			# read in files from there
			allDirFiles <- list.files(targetDir)	

			# files to read in
			toRead <- allDirFiles[grep(".shx", allDirFiles)]

			# the container
			rotatedList <- NULL
			for(i in 1:length(toRead)){
				# read in one bit
				rotatedList[[i]] <- sf::st_read(file.path(targetDir, toRead[i]), quiet=!verbose)
			}

			# make this just one object
			rotated <- NULL
			# the final set of columns
			allColumns <- unique(unlist(lapply(rotatedList, colnames)))
			allColumns <- c(allColumns[allColumns!="geometry"], "geometry")

			for(i in 1:length(rotatedList)){

				# focus on one
				one <- rotatedList[[i]]

				# if columns are missing, add them
				missingColumns <- allColumns[!allColumns%in%colnames(one)]
				if(length(missingColumns)>0){
					for (j in 1:length(missingColumns)){
						one[[missingColumns[j]]] <- NA
					}
				}
				# reorder the columns
				rotated <- rbind(rotated, one[,allColumns])

			}

		}
		
		# if the originals were sf - nothing happens
		if(inherits(originals, "sf") | inherits(originals, "Spatial") | inherits(x, "character")){
			# omission of gplates metadata? 
			if(!gmeta){
				rotated <- 	rotated[, -which(colnames(rotated)%in%c("ANCHOR", "TIME", "FILE1", "RECONFILE1"))]
			}
		}
		
		# but if there was a crs, make sure to change projection to it!
		# is CRS even applicable? 
		if(!is.null(originalCRS)){
			# if it was given in the first place
			if(!is.na(originalCRS)){
				if(verbose) message("Converting back to original CRS.")
				rotated <- sf::st_transform(rotated, originalCRS)
			}
		}
		
	
		# if originals are spatial, make the output spatial
		if(inherits(originals, "Spatial")){
			if(verbose) message("Converting sf to Spatial.")
			# different treatment dependending on whether multiple files or a single one is outpu
			# sp cannot handle heterogeneous geometries
			if(!multiple){
				rotated <- sf::as_Spatial(rotated)
			}else{
				warning("There were multiple reconstruction output, returning as sf.")
			}
		}

		# if originals are not sf and not spatial
		if(!inherits(originals, "sf") & !inherits(originals, "Spatial")){

			# and they are either a matrix or a data frame transform object back to whatever it was
			if((inherits(originals, "matrix") | inherits(originals, "data.frame")) ){
				# run this, except when missing value is indicated
				if(!all(is.na(rotated))){
					# some coordinates probably were missing
					# get the coorindates
					#  reconstructed
					rotatedSF <- rotated
					coords <- sf::st_coordinates(rotated)

					# where to put the new coordinates
					rotated <- originals
					rotated[] <- NA

					# use the IDs to make sure that things really match!
					rownames(rotated) <- theID
					rotated[rotatedSF$ID, ] <- coords

					# copy over the rownames
					rownames(rotated) <- rownames(originals)
				}

			}else{
				if(!is.null(originals)) stop(paste("Unknown output class", class(originals)))
			}
		}


	# 5. Finish
		# remove temporary files
		if(cleanup){		
			if(!inherits(x, "character")){
				system(paste("rm -r ",tempd, "/",layer,"*", sep=""))
			}
		}
	return(rotated)
}

# function to find Gplates installation directory
winDefaultGPlates<-function(){

	# default installation paths
	basic <- c(
		"C:/Program Files/GPlates",
		"C:/Program Files (x86)/GPlates"
		
	)

	versioned <- NULL
	inWhich <- NULL

	# search both possible folders 
	for(i in 1:2){
		# enter program files
		gpver <- list.files(basic[i])
	
		found <- grep("GPlates", gpver)
	
		# grab the latest version
		if(length(found)>0){
			versioned <- gpver[found[length(found)]]
			inWhich <- i
		}
	}
	if(is.null(inWhich)) stop("Could not locate GPlates.")

	# add it 
	dir <- file.path(basic[inWhich], versioned)

	# search executable
	gpfiles <- list.files(dir)

	# grab gplates executable file
	gplat <- grep("gplat",gpfiles)
	potExe <- gpfiles[gplat]
	exe <- potExe[grep("exe",potExe)]

	return(c(dir=dir, exe=exe))
}


macDefaultGPlates <-function(){
# default installation path
	basic <- "/Applications"
	# enter program files
	gpver <- list.files(basic)

	found <- grep("GPlates", gpver)

	# grab the latest version
	if(length(found)>0){
		gpver<- gpver[found[length(found)]]
	}else{
		stop("Could not locate GPlates.")
	}
	dir <-file.path(basic,gpver, "gplates.app/Contents/MacOS")
	exe <-"gplates"
	return(c(dir=dir, exe=exe))
}



testGPlates<- function(gplatesExecutable, verbose){
# ask version
	gplatesTest <- paste(gplatesExecutable, "--v")
	
	# "\"C:/Program Files (x86)/GPlates/GPlates 2.1.0/gplates-2.1.0.exe\" --v"
	# default version
	ver <- NULL
	
	# depending on how much the user wants to see
	if(!verbose){
		opt <- options(show.error.messages = FALSE)
		# revert even if command below fails for some reason
		on.exit(options(opt))

		try(ver <- system(gplatesTest, intern=TRUE,ignore.stdout = TRUE, 
				ignore.stderr = TRUE))
	}else{
		try(ver <- system(gplatesTest, intern=TRUE))
	}
	
	# if gplates is not present
	return(!is.null(ver))
}