Difference between revisions of "OpasnetUtils/Drafts"

library(OpasnetUtils)

ograph <- function( # Määritellään yleisfunktio peruskuvaajan piirtämiseen.
	ovariable, 
	x, 
	y = character(), 
	type = character(), 
	other = character(),
	fill = NA, 
	...
) {
	if(class(ovariable) == "ovariable")  {
		if(nrow(ovariable@output) == 0) ovariable <- EvalOutput(ovariable)
		data <- ovariable@output
		title <- ovariable@name
		if(length(y) == 0) y <- paste(title, "Result", sep = "")
	} else {
		data <- ovariable
		title <- character()
		if(length(y) == 0) y <- "Result"
	}
	if(length(type) == 0) {
		if("Iter" %in% colnames(data)) type <- geom_boxplot() else type <- geom_bar(stat = "identity")
	}
	out <- ggplot(data, aes_string(x = x, y = y, fill = fill)) # Määritellään kuvan sarakkeet
	out <- out + type
	out <- out + theme_grey(base_size=24) # Fontin kokoa suurennetaan
	out <- out + labs(
			title	= title,
			y = paste(unique(data[[paste(title, "Yksikkö", sep = "")]]), sep = "", collapse = ", ")
		)
	out <- out + theme(axis.text.x = element_text(angle = 90, hjust = 1)) # X-akselin tekstit käännetään niin että mahtuvat
	if(length(other) != 0) out <- out + other
	return(out)
}

# fillna takes a data.frame and fills the cells with NA with each level in that column.
# object is the data.frame, marginals is a vector of columns (either column names or positions) that are to be filled.
# This version of fillna accepts column positions (as the previous version) and also column names in marginals.

fillna <- function (object, marginals) {
	a <- dropall(object)
	if(!is.numeric(marginals)) marginals <- match(marginals, colnames(object))
	for (i in marginals) {
		a[[i]] <- as.factor(a[[i]])
		a1 <- a[!is.na(a[[i]]), ]
		a2 <- a[is.na(a[[i]]), ][-i]
		addition <- data.frame(A = levels(a[[i]]))
		colnames(addition) <- colnames(a)[i]
		a2 <- merge(addition, a2)
		a <- rbind(a1, a2)
	}
	return(a)
}

## collapsemarg is a placeholder for a better functionality within CollapseMarginals.
## It takes an ovariable, and summarises all indices in cols using tapply and a user-defined function.
## However, you can also use function "pick" to select locations defined in a list picks found in indices cols.
## Function "unkeep" simply drops the unkept indices without any other operation.
## The output is an ovariable with the same name as the input.
## This was first created for [[:op_fi:Radonin terveysvaikutukset]]

collapsemarg <- function(variable, cols, fun = "sum", picks = list(), ...) { # cols is a character vector, while probs is a list
	out <- dropall(variable@output)
	marginals <- colnames(out)[variable@marginal]

	if(tolower(fun) == "unkeep") { # The function must be a string, otherwise this row will fail.
		out <- out[!colnames(out) %in% cols]
	} else {
		if(tolower(fun) == "pick") {
			for(i in cols) {
				out <- out[out[[i]] %in% picks[[match(i, cols)]] , ]
			}
			cols <- "" # Those locations that were picked are still marginals.
		} else {
			margtemp <- colnames(out)[colnames(out) %in% marginals & !colnames(out) %in% cols] 
			# You must leave at least one index.

			out <- as.data.frame(as.table(tapply(result(variable), out[margtemp], fun)))
			out <- out[!is.na(out$Freq) , ]

			colnames(out)[colnames(out) == "Freq"] <- ifelse(
				length(variable@name) == 0, 
				"Result", 
				paste(variable@name, "Result", sep = "")
			)
		}
	}

	variable@output <- out
	variable@marginal <- colnames(out) %in% marginals & ! colnames(out) %in% cols
	
	return(variable)
}

MyPointKML <- function( # The function creates a KML fille from a SpatialPointsDataFrame object.
	obj = NULL, # Spatial object with the data. A SpatialPointsDataFrame.
	kmlname = "", # Name of the KML fille (does this show on the map?)
	kmldescription = "", # Description of the KML fille (does this show on the map?)
	name = NULL, # Name for each datapoint (vector with the same length as data in obj).
	description = "", # Descrtion of each datapoint (vector with the same length as data in obj).
	icon = "http://maps.google.com/mapfiles/kml/pal4/icon24.png", # Icon shown on pin (?)
	col=NULL # I don't know what this does.
) {

cat("This function MyPointKML is depreciated. Use google.point_kml in OpasnetUtilsExt instead.\n")

    if (is.null(obj)) 
        return(list(header = c("<?xml version=\"1.0\" encoding=\"UTF-8\"?>", 
            "<kml xmlns=\"http://earth.google.com/kml/2.2\">", 
            "<Document>", paste("<name>", kmlname, "</name>", 
                sep = ""), paste("<description><![CDATA[", kmldescription, 
                "]]></description>", sep = "")), footer = c("</Document>", 
            "</kml>")))
    if (class(obj) != "SpatialPointsDataFrame") 
        stop("obj must be of class 'SpatialPointsDataFrame' [package 'sp']")
    if (is.null(name)) {
        name = c()
        for (i in 1:nrow(obj)) name <- append(name, paste("site", 
            i))
    }
    if (length(name) < nrow(obj)) {
        if (length(name) > 1) 
            warning("kmlPoints: length(name) does not match nrow(obj). The first name will be replicated.")
        name <- rep(name, nrow(obj))
    }
    if (length(description) < nrow(obj)) {
        if (length(description) > 1) 
            warning("kmlPoints: length(description) does not match nrow(obj). The first description will be replicated.")
        description <- rep(description, nrow(obj))
    }
    if (length(icon) < nrow(obj)) {
        if (length(icon) > 1) 
            warning("kmlPoints: length(icon) does not match nrow(obj). Only the first one will be used.")
        icon <- icon[1]
    }

	# This is some kind of a colour definition
	
	col2kmlcolor <- function(col) 
		paste(rev(sapply(
			col2rgb(col, TRUE), 
			function(x) sprintf("%02x", x))
		), collapse = "")
	
    kml <- kmlStyle <- ""
    
	# Create the KML fille.
	
	kmlHeader <- c("<?xml version=\"1.0\" encoding=\"UTF-8\"?>","<kml xmlns=\"http://earth.google.com/kml/2.2\">", "<Document>")
    kmlFooter <- c("</Document>", "</kml>")
    
	# Create rows to the KML fille from data in obj.
	
    for (i in 1:nrow(obj)) {
        point <- obj[i, ]
        pt_style <- paste("#style", ifelse(length(icon) == 1, 1, i), sep = "")
        kml <- append(kml, "<Placemark>")
        kml <- append(kml, paste(
		"  <description><![CDATA[",
		name[i],
		": ",
		description[i], 
		"]]></description>", 
		sep = ""
	))
	#kml <- append(kml, "<Style><IconStyle>")
	#kml <- append(kml, paste("<color>", col2kmlcolor(col[i]), "</color>", sep =""))
	#kml <- append(kml, paste("  <Icon><href>", icon, "</href></Icon>", sep = ""))
	#kml <- append(kml, "<scale>0.300000</scale>")
	#kml <- append(kml, "</IconStyle></Style>")
        kml <- append(kml, "  <Point>")
        kml <- append(kml, "    <coordinates>")
        kml <- append(kml, paste(point@coords[1], point@coords[2], sep = ","))
        kml <- append(kml, "    </coordinates>")
        kml <- append(kml, "  </Point>")
        kml <- append(kml, "</Placemark>")
    }
    
    return(paste(paste(c(kmlHeader, kmlStyle, kml, kmlFooter), sep = "", collapse = "\n"), collapse="\n", sep = ""))
}

ova2spat <- function( # This function converts an ovariable into a SpatialPointsDataFrame.
	ovariable, # An evaluated ovariable that has coordinate indices.
	coords, # The names of the coordinate indices as a character vector, first x then y.
	proj4string # Projection identifier or specification as character string. See http://spatialreference.org/
) {
	temp <- ovariable@output
	
	# Transform coordinates into numeric format.
	
	for(i in coords) { 
		if(is(temp[[i]], "factor"))    temp[[i]] <- levels(temp[[i]])[temp[[i]]]
		if(is(temp[[i]], "character")) temp[[i]] <- as.numeric(temp[[i]])
	}
	
	# Define the coordinate points first, then add other ovariable output to it.
	
	sp <- SpatialPoints(temp[coords], CRS(proj4string))
	out <- SpatialPointsDataFrame(sp, temp[!colnames(temp) %in% coords])

	#Transform the projection to longitude-latitude system.
	
	epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
	out <- spTransform(out,epsg4326String)

	return(out)
}

# MyRmap is a function for creating static Google maps as png.
# It is based on MyMap function without the "file destination" parameter
# Requires RgoogleMaps package

MyRmap <- function (
	shp, # a spatial data object
	plotvar, # Name of the column that has the values to be illustrated on the map
	pch = 19, # Shape of the point (19: circle)
	cex = 0.3, # Size of the point
	legend_title = "", # Title of the legend
	legend_position = "topleft",
	numbins = 8, # Number of colour bins in graph
	center, # center of the map
	size = c(640, 480), # size of the map. This produces the right dimensions in Opasnet.
	MINIMUMSIZE = FALSE, 
	RETURNIMAGE = TRUE, 
	GRAYSCALE = FALSE, 
	NEWMAP = TRUE, 
	zoom, 
	verbose = 1, 
	...
) {
	plotvar	<- shp[[plotvar]]
	plotclr	<- brewer.pal(numbins, "Spectral")
	classes	<- classIntervals(plotvar, numbins, style = "quantile")
	colcode	<- findColours(classes, plotclr)
	latR	<- shp@coords[ , 2]
	lonR	<- shp@coords[ , 1]

	#get the bounding box:

	bb <- qbbox(lat = latR, lon = lonR)

	if (missing(zoom)) 
		zoom <- min(MaxZoom(latR, lonR, size))
	if (missing(center)) {
		lat.center <- mean(latR)
		lon.center <- mean(lonR)
	}
	else {
		lat.center <- center[1]
		lon.center <- center[2]
	}
	if (MINIMUMSIZE) {
		ll <- LatLon2XY(latR[1], lonR[1], zoom) # I think the latR and lonR are used here differently than how they
		ur <- LatLon2XY(latR[2], lonR[2], zoom) # are used elsewhere. Thus, if MINIMUMSIZE = TRUE, you may see problems.
		cr <- LatLon2XY(lat.center, lon.center, zoom)
		ll.Rcoords <- Tile2R(ll, cr)
		ur.Rcoords <- Tile2R(ur, cr)
		if (verbose > 1) {
			cat("ll:")
			print(ll)
			print(ll.Rcoords)
			cat("ur:")
			print(ur)
			print(ur.Rcoords)
			cat("cr:")
			print(cr)
		}
		size[1] <- 2 * max(c(ceiling(abs(ll.Rcoords$X)), ceiling(abs(ur.Rcoords$X)))) + 1
		size[2] <- 2 * max(c(ceiling(abs(ll.Rcoords$Y)), ceiling(abs(ur.Rcoords$Y)))) + 1

		if (verbose) cat("new size: ", size, "\n")
	}

	MyMap <- GetMap(
		center = c(lat.center, lon.center),
		zoom = zoom, 
		size = size, 
		RETURNIMAGE = RETURNIMAGE, 
		GRAYSCALE = GRAYSCALE, 
		verbose = verbose, 
		...
	)

	PlotOnStaticMap(MyMap) # Plot an empty map.

	PlotOnStaticMap( # Plot the data points on the map.
		MyMap,
		lat = latR,
		lon = lonR,
		pch = pch,
		cex = cex,
		col = colcode,
		add = T
	)

	legend( # Plot the legend on the map.
		legend_position, 
		legend	= names(attr(colcode, "table")),
		title	= legend_title, 
		fill	= attr(colcode, "palette"), 
		cex	= 1.0, 
		bty	= "y",
		bg	= "white"
	)
}

MyPlotKML <- function(
	shp, # a SpatialPointDataFrame object.
	result = "Result", # The name of  result column in shp.
	rasterization = TRUE, # Whether to rasterize the data or not.
	ncols = 32, # Number or columns in the raster.
	nrows = 32, # Number of rows in the raster.
	fun = mean # function to aggregate data points to the raster.
) {	

cat("Consider merging this function MyPolotKML with google.show_raster_on_maps in OpasnetUtilsExt.\n")

	if(rasterization) {
		#Create blank raster
		rast		<- raster()

		#Set raster extent to that of point data
		extent(rast)	<-extent(shp)

		#Choose number of columns and rows
		ncol(rast)	<- ncols
		nrow(rast)	<- nrows

		#Rasterize point data
		rast2	<- rasterize(shp, rast, shp[[result]], fun = fun)

	}

	start <- 0 # min(shp[[result]])
	end <- max(shp[[result]])
	steps <- approx(c(start,end),n=6)$y
	colors <- rev(rainbow(length(steps), start=0, end=0.50))

	# Create the colorstrip below the map.

	par(mfrow=c(6,1), mar=c(3,1,0,1), cex = 1.5)

	colorstrip <- function(colors, labels)
	{
		count <- length(colors)
		image(
			matrix(1:count, count, 1), 
			col = colors, 
			ylab = "", 
			axes = FALSE
		)
		axis(1,approx(c(0, 1), n=length(labels))$y, labels)
	}

	colorstrip(colors, steps)

	#Plot data

	google.show_raster_on_maps(rast2, col = colors, style = "height:500px;")
}

###### ovaMatrixProduct is a matrix multiplication (%*%) for ovariables.

ovaMatrixProduct <- function(
	e1, # The first ovariable in the multiplication.
	e2, # The second ovariable in the multiplication.
	matrind.a, # The index names used in the matrix (row, column) for the first ovariable. Character vector, length 2.
	matrind.b # The index names used in the matrix (row, column) for the second ovariable. 
) {
	#matrind.a <- c("A", "B")
	#matrind.b <- c("B", "C")

	#Take note of names and marginals
	#Merge the two ovariables
	#Make a data frame with all non-matrix indices, unique rows.
	#Make a for loop and go through non-matrix index rows one by one.
	## Slice the merged ovariable and take the result columns and the matrix indices for each original ovariable.
	## Reshape the two new data.frames to wide format and convert to matrices.
	## Make the matrix product.
	## Convert back to data frame.
	## Append to previous result data.frame together with non-matrix indices.
	# Return a new ovariable.

	matricise <- function( # Produce a matrix of the key data from a data.frame.
		x, # data.frame with the data for making a matrix.
		matrind, # Index names in x for the matrix.
		name # Name of the result column in x.
	) {
		out <- unique(x[c(matrind, name)])

		out <- reshape(out, idvar = matrind[1], timevar = matrind[2], direction = "wide")
		colnames(out) <- gsub(paste(name, ".", sep = ""), "", colnames(out))
		rownames(out) <- out[[1]]
		out[[1]] <- NULL
		out <- as.matrix(out)
		return(out)
	}

	temp <- merge(e1, e2) # Combine ovariables to make all indices match in the future calculations.

	# Create a data.frame that contains unique rows consisting of other marginal indices than those needed in the matrices.

	temp2 <- c(colnames(e1@output)[e1@marginal], colnames(e2@output)[e2@marginal])
	temp2 <- temp2[! temp2 %in% c(matrind.a, matrind.b)]
	nonmatrixind <- unique(temp@output[temp2])

	# Create an output data.frame that collects a matrix product result for each rows of nonmatrixind.

	out <- data.frame()

	for(i in 1:nrow(nonmatrixind))
	{
		temp3 <- merge(temp@output, nonmatrixind[i, , drop = FALSE])

		# Create matrices for ovariables a and b from data.frame temp3.
		
		tempa <- matricise(temp3, matrind.a, comment(result(e1)))	
		tempb <- matricise(temp3, matrind.b, comment(result(e2)))	

		matr <- tempa %*% tempb # The actual matrix product.

		# Convert matrix back to data.frame and combine with other results.
		
		matr <- as.data.frame(as.table(matr))
		colnames(matr) <- c(matrind.a[1], matrind.b[2], "Result")
		matr <- cbind(nonmatrixind[i, , drop = FALSE], matr)
		out <- rbind(out, matr)
	}

	out <- new("ovariable", output = out)
		
	out <- CheckMarginals(out, deps = list(e1, e2), verbose = FALSE)

	return(out)
}

# Merge all but show_bins largest bins of indices cols to 'Other'.

truncateIndex <- function( # Truncates an index to contain only the largest index bins.
	obj, # ovariable to use.
	cols, # names of the columns to truncate.
	bins = 10# Number of bins to show. Other bins will be merged to "Other".
) {
	for(i in 1:length(cols))
	{
		sums <- as.data.frame(as.table(tapply(result(obj), obj@output[cols[i]], sum)))
		bins <- min(bins, length(levels(obj@output[[cols[i]]])))
		limit <- sort(sums$Freq, decreasing = TRUE)[bins]
		keeps <- sums[sums$Freq > limit , cols[i]]
		keeps <- levels(keeps)[keeps]

		levels(obj@output[[cols[i]]]) <- ifelse(
			levels(obj@output[[cols[i]]]) %in% keeps, 
			levels(obj@output[[cols[i]]]), 
			"Other"
		)
	}
	return(obj)
}

findrest <- function (X, cols, total = 1) 
{
	# findrest input is an ovariable that can be integrated over indices cols to result in total.
	# Often it is used with uncertain fractions. The largest fraction is omitted or NA, and it is replaced by whatever is missing from the total.
	if (nrow(X@output) == 0) 
		X <- EvalOutput(X)
	rescol <- comment(result(X))

	# temp is the amount that is still missing from the total.
	temp <- total - oapply(X, cols = cols, FUN = function(x) sum(x, na.rm = TRUE))
	temp@output[[rescol]] <- NULL # Remove old rescol because it would cause trouble later.
	out <- merge(X, temp)@output

	# Replace missing values with values from temp.
	out[[rescol]] <- ifelse(is.na(out[[rescol]]), out$Result, out[[rescol]]) # Result comes from temp
	if(rescol != "Result") out$Result <- NULL
	X@output <- out
	return(X)
}
findrest <- function (X, cols, total = 1)
{
	# findrest input is an ovariable that can be integrated over indices cols to result in total.
	# Often it is used with uncertain fractions. The largest fraction is omitted or NA, and it is replaced by whatever is missing from the total.
	if (nrow(X@output) == 0) 
		X <- EvalOutput(X)
	rescol <- comment(result(X))
	marginals <- colnames(X@output)[X@marginal]

	# temp is the amount that is still missing from the total.
	temp <- total - oapply(X, cols = cols, FUN = function(x) sum(x, na.rm = TRUE))
	temp <- unkeep(temp, cols = rescol) # Remove old rescol because it would cause trouble later.
	colnames(temp@output)[colnames(temp@output) == "Result"] <- "tempResult" 
	temp@name <- "temp" # This is to make sure that merge works.
	print(temp)
	out <- merge(X, temp)@output

	# Replace missing values with values from temp.
	out[[rescol]] <- ifelse(is.na(out[[rescol]]), out$tempResult, out[[rescol]]) # Result comes from temp
	out$tempResult <- NULL
	X@output <- out
	X@marginal <- colnames(X@output) %in% marginals
	return(X)
}

unkeep <- function (X, cols = NA, sources = FALSE, prevresults = FALSE) 
{
	rescol <- paste(X@name, "Result", sep = "")
	if (sources) 
		cols <- c(cols, colnames(X@output)[grepl("Source$", colnames(X@output))])

	for(i in cols) { # Give warning about unkept indices with >1 locations.
		if(
			length(unique(X@output[[i]])) > 1 & 
			X@marginal[match(i, colnames(X@output))]
		) {
			warning(paste(
				"There is >1 unique locations in column", 
				i, ":", 
				paste(as.character(unique(X@output[[i]])), collapse = ", ")
			))
		}
	}

	if (prevresults) 
	cols <- c(cols, colnames(X@output)[grepl("Result$", colnames(X@output)) & 
	colnames(X@output) != rescol])

	marginals <- colnames(X@output)[X@marginal]
	X@output <- X@output[!colnames(X@output) %in% cols]
	X@marginal <- colnames(X@output) %in% marginals

	return(X)
}

objects.store(ograph, fillna, collapsemarg, MyPointKML, ova2spat, MyRmap, MyPlotKML, ovaMatrixProduct, truncateIndex, findrest, unkeep)

cat("The following objects are stored: ograph, fillna, collapsemarg, MyPointKML, ova2spat, MyRmap, MyPlotKML, ovaMatrixProduct, truncateIndex, findrest, unkeep.\n")