Difference between revisions of "OpasnetUtils/Drafts"
From Testiwiki
Jump to: navigation, search
(→Answer) |
(→Answer: MakeOvaDag works with igraph. Some fancy properties) |
||
| Line 28: | Line 28: | ||
MakeOvaDag <- function(...) { | MakeOvaDag <- function(...) { | ||
require(igraph) | require(igraph) | ||
| − | |||
edg <- data.frame() | edg <- data.frame() | ||
| − | for(i in ls(envir = .GlobalEnv)) { | + | ova <- character() |
| + | dec <- character() | ||
| + | deci <- character() | ||
| + | for(i in ls(envir = .GlobalEnv)) { # Find all ovariables | ||
| + | if(class(get(i)) == "ovariable") { | ||
| + | ova <- c(ova, i) | ||
| + | } | ||
| + | if(class(get(i)) == "odecision") { | ||
| + | dec <- c(dec, i) | ||
| + | } | ||
| + | } | ||
| + | for(i in ova) { # Find all edges between ovariables | ||
obj <- get(i) | obj <- get(i) | ||
| − | + | if(nrow(obj@dependencies) > 0) { | |
| − | + | edg <- rbind(edg, data.frame( | |
| − | + | from = obj@dependencies$Name, | |
| − | + | to = obj@name, | |
| − | + | weight = 1 # Not sure yet what should be used as weight. | |
)) | )) | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
} | } | ||
} | } | ||
| − | + | dag <- graph.data.frame(edg, directed = TRUE) # Create directed acyclic graph | |
| − | + | for(i in dec) { | |
| − | # | + | deci <- union(deci, get(i)@dectable$Decision) |
| − | + | } | |
| − | + | for(i in ova) { #Add attributes for vertices | |
| − | dag <- | + | obj <- get(i) |
| − | dag <- dag + vertex( | + | if(!i %in% V(dag)$name) dag <- dag + vertex(obj@name) |
| + | V(dag)$Size[V(dag)$name == i] <- nrow(obj@output) | ||
| + | ind <- colnames(obj@output)[ | ||
| + | obj@marginal & | ||
| + | !grepl("Source$", colnames(obj@output)) & | ||
| + | ! colnames(obj@output) %in% deci | ||
| + | ] | ||
| + | decii <- deci[deci %in% colnames(obj@output)] | ||
| + | if(length(ind) > 0) V(dag)$Indices[V(dag)$name == i] <- paste(ind, collapse = ", ") | ||
| + | if(length(decii) > 0) V(dag)$Decisions[V(dag)$name == i] <- paste(decii, collapse = ", ") | ||
| + | } | ||
| + | V(dag)$Size[is.na(V(dag)$Size)] <- 1 | ||
| + | V(dag)$Size[V(dag)$Size == 0] <- 1 | ||
| + | V(dag)$Indices[is.na(V(dag)$Indices)] <- "" | ||
| + | V(dag)$Decisions[is.na(V(dag)$Decisions)] <- "" | ||
| + | |||
| + | plot( | ||
| + | dag, | ||
| + | vertex.label.cex = 0.8, | ||
| + | vertex.size = log(V(dag)$Size)+2, # Vertex size is (non-linearly) relative to rows in output. | ||
| + | vertex.color = ifelse(nchar(V(dag)$Decisions) > 0, "Red", "SkyBlue2"), | ||
| + | layout = layout.fruchterman.reingold | ||
| + | ) | ||
return(dag) | return(dag) | ||
} | } | ||
Revision as of 14:21, 21 September 2015
This page is a method.
The page identifier is Op_en6007 |
|---|
| Moderator:Jouni (see all) |
|
Give your opinion to the peer rating of the content of this page. |
| Upload data
|
Question
Which functions are so useful that they should be taken into OpasnetUtils package? This page contains draft function which will be included when they are good enough and found important.
Answer
Call the objects stored by this code from another rode with this command:
objects.latest("Op_en6007", code_name = "answer")
library(OpasnetUtils)
##### Makes a directed acyclic graph (DAG) out of all ovariables and their dependencies in the global environment.
##### Requires igraph package.
### ideas of additional atteibutes: decision indices are shown as specific colour on vertex border.
## voi analysis results are shown as the size of vertex. Can this be retrieved from somewhere automatically? Only if there is a specific object eg voi.
## good functions for layout: layout.fruchterman.reingold() layout.reingold.tilford() layout.lgl()
MakeOvaDag <- function(...) {
require(igraph)
edg <- data.frame()
ova <- character()
dec <- character()
deci <- character()
for(i in ls(envir = .GlobalEnv)) { # Find all ovariables
if(class(get(i)) == "ovariable") {
ova <- c(ova, i)
}
if(class(get(i)) == "odecision") {
dec <- c(dec, i)
}
}
for(i in ova) { # Find all edges between ovariables
obj <- get(i)
if(nrow(obj@dependencies) > 0) {
edg <- rbind(edg, data.frame(
from = obj@dependencies$Name,
to = obj@name,
weight = 1 # Not sure yet what should be used as weight.
))
}
}
dag <- graph.data.frame(edg, directed = TRUE) # Create directed acyclic graph
for(i in dec) {
deci <- union(deci, get(i)@dectable$Decision)
}
for(i in ova) { #Add attributes for vertices
obj <- get(i)
if(!i %in% V(dag)$name) dag <- dag + vertex(obj@name)
V(dag)$Size[V(dag)$name == i] <- nrow(obj@output)
ind <- colnames(obj@output)[
obj@marginal &
!grepl("Source$", colnames(obj@output)) &
! colnames(obj@output) %in% deci
]
decii <- deci[deci %in% colnames(obj@output)]
if(length(ind) > 0) V(dag)$Indices[V(dag)$name == i] <- paste(ind, collapse = ", ")
if(length(decii) > 0) V(dag)$Decisions[V(dag)$name == i] <- paste(decii, collapse = ", ")
}
V(dag)$Size[is.na(V(dag)$Size)] <- 1
V(dag)$Size[V(dag)$Size == 0] <- 1
V(dag)$Indices[is.na(V(dag)$Indices)] <- ""
V(dag)$Decisions[is.na(V(dag)$Decisions)] <- ""
plot(
dag,
vertex.label.cex = 0.8,
vertex.size = log(V(dag)$Size)+2, # Vertex size is (non-linearly) relative to rows in output.
vertex.color = ifelse(nchar(V(dag)$Decisions) > 0, "Red", "SkyBlue2"),
layout = layout.fruchterman.reingold
)
return(dag)
}
#dag <- MakeOvaDag()
#plot(dag)
##################### Forgets decisions so that decision indices will be recreated.
forgetDecisions <- function() {
for(i in ls(envir = openv)) {
if("dec_check" %in% names(openv[[i]])) openv[[i]]$dec_check <- FALSE
}
return(cat("Decisions were forgotten.\n"))
}
#################### Perus-ggplot ovariableille
oggplot <- function (
ova, # ovariable to be plotted
x, # Index for x axis
weight = NULL, # Index for y axis (default: result column)
fill = NULL, # Index for colour code
base_size = BS, # Base size for graph font (object BS must exist!)
turnx = FALSE, # Turn x axis labels vertically?
binwidth = NULL # Width of bins
)
{
if(exists("suomenna")) ova <- suomenna(ova)
if (is.null(weight))
weight <- paste(ova@name, "Result", sep = "")
if("Iter" %in% colnames(ova@output))
ova <- oapply(ova, cols = "Iter", FUN = mean)
plo <- ggplot(ova@output, aes_string(x = x, fill = fill, weight = weight)) +
theme_gray(base_size = base_size)
if(any(ova@output[[weight]] >= 0)) {
plo <- plo + geom_bar(
data = subset(ova@output, ova@output[[weight]] >= 0),
position = "stack",
binwidth = binwidth
)
}
if(any(ova@output[[weight]] < 0)) {
plo <- plo + geom_bar(
data = subset(ova@output, ova@output[[weight]] < 0),
position = "stack",
binwidth = binwidth
) + geom_hline(aes(yintercept = 0))
}
if (turnx)
plo <- plo + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
return(plo)
}
################## Sähkön hinta tunneittain
#price <- opbase.data(ident="op_en7353")
#temperature <- opbase.data("op_en6315.2014_5_2015")
#temperature$Date <- substr(temperature$Date, 0, 11)
#price$Date <- substr(price$Date, 0, 11)
#mon <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
#for (i in mon) {
# price$Date <- gsub(i, which(mon == i), as.character(price$Date))
#}
#for (i in mon) {
# temperature$Date <- gsub(i, which(mon == i), as.character(temperature$Date))
#}
#price$Hours <- substr(price$Hours, 0, 2)
#price$Hours <- paste(price$Hours, ":00:00", sep="")
#temperature$Time <- paste(temperature$Time, ":00", sep="")
#as.character(temperature$Result)
#as.numeric(temperature$Result)
#cut(temperature$Result, breaks = c(-21, -18, -15, -12, -9, -6, -3, 0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30), #include.lowest=TRUE)
#DateTime <- as.POSIXct(paste(temperature$Date, temperature$Time), format="%Y-%m-%d %H:%M:%S")
#DateHours <- as.POSIXct(paste(price$Date, price$Hours), format="%Y-%m-%d %H:%M:%S")
################## Suomentaja
suomenna <- function(ova) {
if(exists("finnish")) if(finnish) {
d <- ova@output
colnames(d) <- gsub("[ \\.]", "_", colnames(d))
if("Decision_maker" %in% colnames(d)) {
d$Decision_maker <- as.factor(d$Decision_maker)
levels(d$Decision_maker)[levels(d$Decision_maker) == "Builders"] <- "Rakennuttaja"
levels(d$Decision_maker)[levels(d$Decision_maker) == "Building owner"] <- "Rakennuksen omistaja"
}
if("Decision" %in% colnames(ova@output)) {
levels(d$Decision)[levels(d$Decision) == "EnergySavingPolicy"] <- "Energiansäästöpolitiikka"
levels(d$Decision)[levels(d$Decision) == "PlantPolicy"] <- "Voimalapolitiikka"
}
if("Option" %in% colnames(ova@output)) {
levels(d$Option)[levels(d$Option) == "BAU"] <- "Tätä menoa"
levels(d$Option)[levels(d$Option) == "Energy saving moderate"] <- "Kohtuullinen energiansäästö"
levels(d$Option)[levels(d$Option) == "Energy saving total"] <- "Täysi energiansäästö"
}
if("Building" %in% colnames(ova@output)) {
levels(d$Building)[levels(d$Building) == "Apartment houses"] <- "Kerrostalot"
levels(d$Building)[levels(d$Building) == "Commercial"] <- "Kaupalliset"
levels(d$Building)[levels(d$Building) == "Detached houses"] <- "Omakotitalot"
levels(d$Building)[levels(d$Building) == "Educational"] <- "Opetusala"
levels(d$Building)[levels(d$Building) == "Health and social sector"] <- "Terveys- ja sosiaaliala"
levels(d$Building)[levels(d$Building) == "Industrial"] <- "Teollisuus"
levels(d$Building)[levels(d$Building) == "Leisure houses"] <- "Mökki"
levels(d$Building)[levels(d$Building) == "Offices"] <- "Toimistot"
levels(d$Building)[levels(d$Building) == "Other"] <- "Muu"
levels(d$Building)[levels(d$Building) == "Public"] <- "Julkinen"
levels(d$Building)[levels(d$Building) == "Row houses"] <- "Rivitalot"
levels(d$Building)[levels(d$Building) == "Sports"] <- "Urheilu"
}
if("Efficiency" %in% colnames(ova@output)) {
levels(d$Efficiency)[levels(d$Efficiency) == "Traditional"] <- "Perinteinen"
levels(d$Efficiency)[levels(d$Efficiency) == "Old"] <- "Vanha"
levels(d$Efficiency)[levels(d$Efficiency) == "New"] <- "Uusi"
levels(d$Efficiency)[levels(d$Efficiency) == "Low-energy"] <- "Matalaenerginen"
levels(d$Efficiency)[levels(d$Efficiency) == "Passive"] <- "Passiivitalo"
}
if("Renovation" %in% colnames(ova@output)) {
levels(d$Renovation)[levels(d$Renovation) == "None"] <- "Ei mitään"
levels(d$Renovation)[levels(d$Renovation) == "General"] <- "Yleinen"
levels(d$Renovation)[levels(d$Renovation) == "Windows"] <- "Ikkunat"
levels(d$Renovation)[levels(d$Renovation) == "Techical systems"] <- "Tekniset"
levels(d$Renovation)[levels(d$Renovation) == "Sheath reform"] <- "Seinät ja katto"
}
if("Plant" %in% colnames(ova@output)) {
levels(d$Plant)[levels(d$Plant) == "Biofuel heat plants"] <- "Biolämpölaitokset"
levels(d$Plant)[levels(d$Plant) == "CHP diesel generators"] <- "CHP dieselgeneraattorit"
levels(d$Plant)[levels(d$Plant) == "Deep-drill heat"] <- "Syväporattu maalämpö"
levels(d$Plant)[levels(d$Plant) == "Domestic"] <- "Omaan käyttöön"
levels(d$Plant)[levels(d$Plant) == "Hanasaari"] <- "Hanasaari"
levels(d$Plant)[levels(d$Plant) == "Hanasaari biofuel renovation"] <- "Hanasaari bio"
levels(d$Plant)[levels(d$Plant) == "Household air heat pumps"] <- "Kotitalouden ilmalämpöpumppu"
levels(d$Plant)[levels(d$Plant) == "Household air conditioning"] <- "Kotitalouden ilmastointi"
levels(d$Plant)[levels(d$Plant) == "Household geothermal heat"] <- "Kotitalouden maalämpö"
levels(d$Plant)[levels(d$Plant) == "Household solar"] <- "Kotitalouden aurinkovoima"
levels(d$Plant)[levels(d$Plant) == "Katri Vala cooling"] <- "Katri Vala viilennys"
levels(d$Plant)[levels(d$Plant) == "Katri Vala heat"] <- "Katri Vala lämpö"
levels(d$Plant)[levels(d$Plant) == "Kellosaari back-up plant"] <- "Kellosaaren varavoimala"
levels(d$Plant)[levels(d$Plant) == "Kymijoki River's plants"] <- "Kymijoen vesivoimalat"
levels(d$Plant)[levels(d$Plant) == "Loviisa nuclear heat"] <- "Loviisan ydinvoimalämpö"
levels(d$Plant)[levels(d$Plant) == "Neste oil refinery heat"] <- "Nesteen öljyjalostamolämpö"
levels(d$Plant)[levels(d$Plant) == "Other"] <- "Muu"
levels(d$Plant)[levels(d$Plant) == "Salmisaari A&B"] <- "Salmisaari A&B"
levels(d$Plant)[levels(d$Plant) == "Salmisaari biofuel renovation"] <- "Salmisaari biokorjaus"
levels(d$Plant)[levels(d$Plant) == "Sea heat pump"] <- "Merilämpö"
levels(d$Plant)[levels(d$Plant) == "Sea heat pump for cooling"] <- "Merilämpö viilennykseen"
levels(d$Plant)[levels(d$Plant) == "Small-scale wood burning"] <- "Puun pienpoltto"
levels(d$Plant)[levels(d$Plant) == "Small fuel oil heat plants"] <- "Pienet öljylämpölaitokset"
levels(d$Plant)[levels(d$Plant) == "Small gas heat plants"] <- "Pienet kaasulämpölaitokset"
levels(d$Plant)[levels(d$Plant) == "Suvilahti power storage"] <- "Suvilahden voimavaraaja"
levels(d$Plant)[levels(d$Plant) == "Suvilahti solar"] <- "Suvilahden aurinkovoima"
levels(d$Plant)[levels(d$Plant) == "Vuosaari A&B"] <- "Vuosaari A&B"
levels(d$Plant)[levels(d$Plant) == "Vuosaari C biofuel"] <- "Vuosaari C bio"
levels(d$Plant)[levels(d$Plant) == "Wind mills"] <- "Tuulivoimalat"
}
if("Fuel" %in% colnames(ova@output)) {
levels(d$Fuel)[levels(d$Fuel) == "Electricity"] <- "Sähkö"
levels(d$Fuel)[levels(d$Fuel) == "Heat"] <- "Lämpö"
levels(d$Fuel)[levels(d$Fuel) == "Biofuel"] <- "Biopolttoaine"
levels(d$Fuel)[levels(d$Fuel) == "Coal"] <- "Kivihiili"
levels(d$Fuel)[levels(d$Fuel) == "Fuel oil"] <- "Polttoöljy"
levels(d$Fuel)[levels(d$Fuel) == "Gas"] <- "Maakaasu"
levels(d$Fuel)[levels(d$Fuel) == "Light oil"] <- "Kevytöljy"
levels(d$Fuel)[levels(d$Fuel) == "Wood"] <- "Puu"
}
if("Heating" %in% colnames(ova@output)) {
levels(d$Heating)[levels(d$Heating) == "District"] <- "Kaukolämpö"
levels(d$Heating)[levels(d$Heating) == "Electricity"] <- "Sähkölämmitys"
levels(d$Heating)[levels(d$Heating) == "Oil"] <- "Öljy"
levels(d$Heating)[levels(d$Heating) == "Other"] <- "Muu"
}
colnames(d)[colnames(d) == "Decision_maker"] <- "Päätöksentekijä"
colnames(d)[colnames(d) == "Decision"] <- "Päätös"
colnames(d)[colnames(d) == "Option"] <- "Vaihtoehto"
colnames(d)[colnames(d) == "Building"] <- "Rakennus"
colnames(d)[colnames(d) == "Efficiency"] <- "Tehokkuus"
colnames(d)[colnames(d) == "Renovation"] <- "Korjaukset"
colnames(d)[colnames(d) == "Plant"] <- "Voimala"
colnames(d)[colnames(d) == "Fuel"] <- "Polttoaine"
ova@output <- d
}
return(ova)
}
ograph <- function( # Määritellään yleisfunktio peruskuvaajan piirtämiseen.
ovariable,
x,
y = character(),
type = character(),
other = character(),
fill = NA,
...
) {
cat("This function ograph is depreciated. Use oggplot instead.\n")
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.
# fillna was updated in OpasnetUtils and therefore removed from here.
## 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;")
}
# Merge all but show_bins largest bins of indices cols to 'Other'.
truncateIndex <- function( # Truncates indices to contain only the largest index bins.
obj, # ovariable to use.
cols, # names of the columns to truncate.
bins = rep(10, length(cols)), # Number of bins to show, including Others. Smallest locations will be lumped to bin "Other".
sum_others = TRUE # Should "Other" be summed to maintain marginal status
) {
if(nrow(obj@output) == 0) stop("Ovariable ", obj@name, " not evaluated.\n")
test <- oapply(abs(obj), INDEX = cols, sum, na.rm = TRUE)
if(length(cols) > 1 & length(bins) == 1) bins <- rep(bins, length(cols))
for(i in 1:length(cols))
{
test2 <- oapply(test, INDEX = cols[i], sum)
test2@output <- test2@output[result(test2) > 0 , ]
temp <- as.factor(obj@output[[cols[i]]])
location_weight_order <- order(result(test2), decreasing = TRUE)
keeps <- test2@output[[cols[i]]][location_weight_order[0:min(bins[i] - 1, nrow(test2@output))]]
levels(temp)[!levels(temp) %in% keeps] <- "Other"
temp <- factor(temp, levels = c(levels(temp)[levels(temp) != "Other"], "Other"))
obj@output[[cols[i]]] <- temp
}
# After changing some locations to "Other", sum along indices to avoid problems
if(sum_others) {
obj <- oapply(obj, cols = "", FUN = sum, na.rm = TRUE)
}
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. One (often 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 <- paste(X@name, "Result", sep = "")
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))
# Remove old rescol because it would cause trouble later.
if(rescol != paste(temp@name, "Result", sep = "")) temp <- unkeep(temp, cols = rescol)
colnames(temp@output)[colnames(temp@output) == paste(temp@name, "Result", sep = "")] <- "tempResult"
temp@name <- "temp" # This is to make sure that merge works.
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)
}
timing <- function(dat, timecol = NA, weeks = 6, tz = "EET")
# timing converts character or numeric inputs into times using a few default formats
{
if(is.data.frame(dat)) # Turn dat into a data.frame in all cases.
{
timesall <- (dat[timecol])
} else {
timesall <- data.frame(Timecol = dat)
timecol <- "Timecol"
}
temprow <- character() # A whole row of timesall collapsed into a string.
for(i in 1:nrow(timesall))
{
temprow[i] <- tolower(paste(t(timesall)[ , i], collapse = ""))
}
weekdays <- data.frame(
Name = c("su", "ma", "ti", "ke", "to", "pe", "la",
"sun", "mon", "tue", "wed", "thu", "fri", "sat",
"sön", "mon", "tis", "ons", "tor", "fre", "lör"),
Number = rep(0:6, times = 3)
)
# day is the number of the weekday in case of repeating events.
day <- rep(NA, nrow(timesall))
for(i in 1:nrow(weekdays))
{
temp1 <- grepl(weekdays$Name[i], temprow) # Is any weekday mentioned in temprow?
day <- ifelse(is.na(day) & temp1, weekdays$Number[i], day) # Find the weekday number.
}
# repall are repeating events, x are non-repeating events.
repall <- timesall[!is.na(day) , ]
xall <- timesall[is.na(day) , ]
starting <- NA
########### First, change non-repeating timedates.
# Change timedate into POSIXct assuming formats 15.3.2013 or 2013-03-15 and 15:24 or 15.24.
# Note! 13 and 2013 mean 1.1.2013 and 3.2013 means 1.3.2013
if(nrow(xall) > 0)
{
xout <- data.frame(Datrow = (1:nrow(timesall))[is.na(day)]) # Row numbers from dat.
for(j in colnames(xall))
{
x <- xall[[j]]
if(is.factor(x)) x <- levels(x)[x]
x <- ifelse(grepl("^[0-9][0-9]$", x), paste("20", x, sep = ""), x)
x <- ifelse(grepl("^[0-9][0-9][0-9][0-9]$", x), paste("1.1.", x, sep = ""), x)
x <- ifelse(grepl("^[0-9].[0-9][0-9][0-9][0-9]$", x), paste("1.", x, sep = ""), x)
x <- ifelse(grepl("^[0-9][0-9].[0-9][0-9][0-9][0-9]$", x), paste("1.", x, sep = ""), x)
temp <- x
x <- as.POSIXct(temp, format = "%d.%m.%Y %H:%M", tz = tz)
x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%d.%m.%Y %H.%M", tz = tz))
x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%Y-%m-%d %H:%M", tz = tz))
x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%Y-%m-%d %H.%M", tz = tz))
x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%d.%m.%Y", tz = tz))
x <- ifelse(!is.na(x), x, as.POSIXct(temp, format = "%d.%m.%y", tz = tz))
xout <- cbind(xout, X = as.POSIXct(x, origin = "1970-01-01", tz = tz))
# Intermediate values turn into numeric, therefore turned back to POSIX.
starting <- min(c(starting, xout$X), na.rm = TRUE)
colnames(xout)[colnames(xout) == "X"] <- j
}
}
############## Then, change weekly repeating timedates.
# repall must have format ma 9:00 or TUESDAY 8.24. Weekday is case-insensitive and can be abbrevieated.
# The start and end times are assumed to be on the same day. The name of day can be on either column.
if(nrow(repall) > 0)
{
for(j in colnames(repall))
{
reptime <- gsub("[[:alpha:] ]", "", repall[[j]]) # Remove alphabets and spaces.
reptime <- gsub("\\.", ":", reptime)
reptime <- strsplit(reptime, split = ":")
temp2 <- numeric()
for(i in 1:length(reptime))
{
temp2[i] <- as.numeric(reptime[[i]][1]) * 3600 + as.numeric(reptime[[i]][2]) * 60
}
reptime <- temp2
if(is.na(starting)) starting <- paste(format(Sys.Date(), format = "%Y"), "-01-01", sep = "")
starting <- as.POSIXlt(starting, origin = "1970-01-01") # First day of year.
starting$mday <- starting$mday - as.numeric(format(starting, format = "%w")) + day[!is.na(day)] # Previous Sunday plus weekdaynumber.
reps <- data.frame()
temp3 <- starting
for(i in 1:weeks)
{
reps <- rbind(reps, data.frame(
Datrow = (1:nrow(timesall))[!is.na(day)],
X = as.POSIXct(temp3, origin = "1970-01-01", tz = tz) + reptime
))
temp3$mday <- temp3$mday + 7 # Make a weekly event.
}
colnames(reps)[colnames(reps) == "X"] <- j
if(j == colnames(repall)[1]) repout <- reps else repout <- cbind(repout, reps[j])
}
}
out <- rbind(xout, repout)
if(is.data.frame(dat))
{
dat$Datrow <- 1:nrow(dat)
dat <- dat[!colnames(dat) %in% timecol]
out <- merge(dat, out)
out <- out[colnames(out) != "Datrow"]
} else {
out <- out[[timecol]]
}
return(out)
}
# Funktio makeTimeline ottaa tapahtumalistauksen ja rakentaa siitä aikajanan. Parametrit:
# event = data.framena tapahtumalistaus, joka sisältää ainakin alku- ja loppuajan (Alku, Loppu) ja
# tapahtumatiedon sekä mahdollisesti toiston ja keston (Toisto = aikaväli päivinä,
# Kesto = viimeinen tapahtuma-aika).
# timeformat = jos TRUE, oletetaan POSIX-muotoiseksi ja muutetaan operointia varten sekunneiksi.
# Jos FALSE, oletetaan reaaliluvuksi ja operoidaan suoraan luvuilla.
makeTimeline <- function(event, timeformat = TRUE) {
# eventiin luodaan Toisto ja Kesto, jos jompikumpi puuttuu.
if(!all(c("Toisto", "Kesto") %in% colnames(event)))
{
event$Toisto <- NA
event$Kesto <- NA
}
if(timeformat)
{
for(m in c("Alku", "Loppu", "Kesto"))
{
event[[m]] <- as.double(as.POSIXct(event[[m]])) # Muutetaan aika sekunneiksi.
}
event$Toisto <- event$Toisto * 3600 * 24 # Muutetaan Toisto päivistä sekunneiksi.
}
# Jos on puuttuvaa tietoa kestosta tai toistosta, korvataan inertillä datalla.
test <- event$Toisto == 0 | event$Kesto == 0 | is.na(event$Toisto) | is.na(event$Kesto)
event$Toisto <- ifelse(test, 1, event$Toisto)
event$Kesto <- ifelse(test, event$Alku, event$Kesto)
# Luodaan aikajanan alkupiste. Eventrow otetaan aikanaan tyhjältä riviltä.
timeline <- data.frame(Time = min(event$Alku), Eventrow = nrow(event)+1)
for(i in 1:nrow(event)) { # Toista jokaiselle havaintoriville
times <- 0:floor((event$Kesto[i] - event$Alku[i]) / event$Toisto[i]) # Toistojen määrä
# Toista tapahtumaa Kestoon asti.
temp <- data.frame(
Time = event$Alku[i] + times * event$Toisto[i],
End = event$Loppu[i] + times * event$Toisto[i],
EventrowStart = i
)
timeline <- merge(timeline, temp[ , c("Time", "EventrowStart")], all = TRUE) # Lisätään tapahtumat aikajanaan.
colnames(temp) <- c("Remove", "Time", "EventrowEnd") # Muutetaan otsikot, koska nyt halutaan mergata loppuhetket aikajanaan.
timeline <- merge(timeline, temp[, c("Time", "EventrowEnd")], all = TRUE)
for(j in 2:nrow(timeline)) {
# Tämä luuppi käy aikajanan läpi ja täydentää tapahtumat.
# Ensin kaikkiin uusiin aikapisteisiin jatketaan sitä aiempaa toimintaa, joka oli menossa edellisessä pisteessä.
if(is.na(timeline$Eventrow[j]))
{
timeline$Eventrow[j] <- timeline$Eventrow[j-1]
}
# Sitten jatketaan uutta toimintaa niihin uusiin pisteisiin, jotka eivät ole loppupisteitä.
if(is.na(timeline$EventrowStart[j]) & is.na(timeline$EventrowEnd[j]))
{
timeline$EventrowStart[j] <- timeline$EventrowStart[j-1]
}
}
# Jos uutta toimintaa on olemassa, statuksena käytetään sitä, muutoin statusta eli aiempaa toimintaa.
timeline$Eventrow <- ifelse(!is.na(timeline$EventrowStart), timeline$EventrowStart, timeline$Eventrow)
# Leikataan turhat sarakkeet pois ja siirrytään seuraavalle event-riville
timeline <- timeline[ , c("Time", "Eventrow")]
}
event <- rbind(event, rep(NA, ncol(event))) # Lisää rivi loppuhetkeä varten
event$Alku[nrow(event)] <- max(timeline$Time)
event$Eventrow <- row(event)[ , 1]
timeline <- merge(timeline, event) # Yhdistetään Statuksen eli eventin rivinumeron avulla.
timeline <- timeline[!colnames(timeline) %in% c("Eventrow", "Alku", "Loppu", "Toisto", "Kesto")]
timeline <- timeline[order(timeline$Time) , ]
if(timeformat) timeline$Time <- as.POSIXct(timeline$Time, origin = "1970-01-01")
return(timeline)
}
# Calculate the cumulative impact of the events on building stock to given years
timepoints <- function(
# Function timepoints takes an event list and turns that into existing crosscutting situations at
# timepoints defined by years. The output will have index Time.
# In other words, this will integrate over obstime at specified timepoints.
X, # X must be an ovariable with a column of the same name as obstime.
obstime, # obstime must be a single-column data.frame of observation times.
sumtimecol = TRUE # Should the timecol be summed up?
# obstime and timecol may be numeric (by coercion) or POSIXt.
) {
timecol <- colnames(obstime)
marginals <- colnames(X@output)[X@marginal]
# tapply (and therefore possibly oapply) changes continuous indices to factors! Must change back by hand.
if("factor" %in% c(class(obstime[[timecol]]), class(X@output[[timecol]]))) {
X@output[[timecol]] <- as.numeric(as.character(X@output[[timecol]]))
obstime[[timecol]] <- as.numeric(as.character(obstime[[timecol]]))
}
out <- data.frame()
if(sumtimecol) by <- setdiff(marginals, timecol) else by <- marginals
for(i in obstime[[timecol]]) {
temp <- X@output[X@output[[timecol]] <= i , ]
if(nrow(temp) > 0) {
temp <- aggregate(
temp[paste(X@name, "Result", sep = "")],
by = temp[colnames(temp) %in% by],
FUN = sum
)
}
if(nrow(temp) > 0) out <- rbind(out, data.frame(Time = i, temp))
}
X@output <- out
X@marginal <- colnames(out) %in% c("Time", marginals) # Add Time to marginal
return(X)
}
# ana2ova takes in variable tables from Analytica (produced with Copy Table). The output is a data.frame in long format.
ana2ova <- function(dat)
{
i <- 1 # Number of indices
cols <- character() # Names of indices
locs <- list() # Vectors of locations for each index.
out <- data.frame() # output
repeat{ # Find out cols and locs from dat.
temp <- strsplit(dat[i + 1], split = "\t")[[1]]
cols[i] <- temp[1]
if(length(temp) == 1) break
locs[[i]] <- temp[2:length(temp)]
i <- i + 1
}
locs[[i]] <- NA # The innermost index should have its place in locs even if the locations are not yet known.
inner <- (1:length(dat))[dat %in% cols[i]] # Places where the innermost index is mentioned (the data starts from the next row).
if(!all(dat[inner - 1] == dat[i]) & i > 1) stop("Structure incorrect\n")
if(length(inner) == 1) # len is the length of the innermost index.
{
len <- length(dat) - i - 1
} else {
len <- inner[2] - inner[1] - i
}
for(k in inner) # Go through each table with the inner and second most inner indices.
{
ins <- strsplit(dat[(k + 1):(k + len)], split = "\t") # Make a data.frame.
ins <- data.frame(matrix(unlist(ins), nrow = len, byrow = TRUE))
if(i == 1) cn <- "value" else cn <- locs[[i - 1]]
colnames(ins) <- c(cols[i], cn) # Give location names to columns
if(i > 1) ins <- melt(ins, id.vars = cols[i], variable.name = cols[i - 1]) # If several columns.
if(i > 2) # If there are several tables.
{
uplocs <- data.frame(Removethis = 0)
for(l in 1:(i - 2))
{
uplocs[[cols[l]]] <- strsplit(dat[k - i + l], split = "\t")[[1]][2]
}
out <- rbind(out, cbind(uplocs, ins)) # Collects all tables into a data.frame.
} else {
out <- ins
}
}
out$Removethis <- NULL
colnames(out)[colnames(out) == "value"] <- "Result"
return(out)
}
orbind2 <- function( # Like orbind but the value is an ovariable.
o1, # ovariable whose slots are used in the value.
o2, # ovariable
use_fillna = FALSE, # Do we use fillna to fill in the NA values in indices?
warn = "" # What warning is given if fillna is used?
) {
x <- unkeep(o1 * 1, prevresults = TRUE, sources = TRUE)
y <- unkeep(o2 * 1, prevresults = TRUE, sources = TRUE)
xmarg <- colnames(x@output)[x@marginal]
ymarg <- colnames(y@output)[y@marginal]
for(i in xmarg) x@output[[i]] <- as.factor(x@output[[i]])
for(i in ymarg) y@output[[i]] <- as.factor(y@output[[i]])
out <- o1
out@output <- orbind(x, y)
if(use_fillna) {
b <- character()
for(i in colnames(out@output)[out@marginal]) {if(any(is.na(out@output[[i]]))) b <- c(b, i)}
if(length(b) > 0) {
out@output <- fillna(out@output, b)
warning(warn, "\nMissing values had to be filled by function fillna in indices: ", b, "\n")
}
}
colnames(out@output)[colnames(out@output) == "Result"] <- paste(o1@name, "Result", sep = "")
out@marginal <- colnames(out@output) %in% c(xmarg, ymarg)
return(out)
}
objects.store(list = ls())
cat(paste("All objects in the global namespace were stored.\n"))
|
See also
- OpasnetUtils/Ograph, a previous code, now depreciated.
- en:Matrix multiplication in Wikipedia, Matmult in R
References
Related files
<mfanonymousfilelist></mfanonymousfilelist>
