Difference between revisions of "OpasnetUtils/Drafts"

From Testiwiki
Jump to: navigation, search
(→‎Answer: ova2spat updated)
(→‎Calculations: showind works now)
 
(8 intermediate revisions by the same user not shown)
Line 13: Line 13:
 
Call the objects stored by this code from another rode with this command:
 
Call the objects stored by this code from another rode with this command:
  
−
  objects.latest("Op_en6007", code_name = "answer")
+
  objects.latest("Op_en6007", code_name = "answer") # Old version that fetches all objects, depreciated and not updated.
 +
objects.latest("Op_en6007", code_name = "diagnostics") # Functions for ovariable and model diagnostics: ovashapetest, showLoctable, binoptest
 +
objects.latest("Op_en6007", code_name = "webropol") # Functions for operating with Webropol data
 +
objects.latest("Op_en6007", code_name = "miscellaneous") # Functions for various tasks
 +
objects.latest("Op_en6007", code_name = "gis") # Functions for ovariable, KML and Googl maps interactions
  
−
<rcode name="answer" embed=1 store=1>
+
== Rationale ==
 +
 
 +
=== Calculations ===
 +
 
 +
'''Functions for ovariable diagnostics
 +
showind has problems with get() but this version of code was acceptable [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=QcfKMkCd2ewUtZqP].
 +
 
 +
<rcode name="diagnostics" embed=0>
 +
#This is code Op_en6007/diagnostics on page [OpasnetUtils/Drafts]]
 +
 
 +
library(OpasnetUtils)
 +
 
 +
# Shows a table about ovariables and their index and location changes compared with parents.
 +
# showind has problems with get().
 +
showind <- function(name = ".GlobalEnv", sources = FALSE, prevresults = FALSE) {
 +
  # i ovariable
 +
  # k parent ovariable
 +
  # l index in (parent) ovariable
 +
  deptable <- data.frame()
 +
  for(i in ls(name = name)) {
 +
    d = list(get(i))[[1]]
 +
    if(class(d) == "ovariable") {
 +
      depind <- list()
 +
      if(nrow(d@dependencies)>0) {
 +
        dep <- paste(d@dependencies$Name, collapse = ", ")
 +
        for(k in d@dependencies$Name){
 +
          if(!exists(k)) cat(k, "does not exist.\n") else {
 +
            if(class(get(k)) != "ovariable") cat(k, "is not an ovariable.\n") else {
 +
              ko <- list(get(k)@output)[[1]]
 +
              if("Iter" %in% colnames(ko)) ko$Iter <- as.factor(max(as.numeric(as.character(ko$Iter))))
 +
              cols <- colnames(ko)
 +
              if(!sources) cols <- cols[!grepl("Source$", cols)]
 +
              if(!prevresults) cols <- cols[!grepl("Result$", cols)]
 +
              for(l in cols) {
 +
                if(l %in% names(depind)) {
 +
                  depind[[l]] <- union(depind[[l]], unique(ko[[l]]))
 +
                } else {
 +
                  newind <- list(unique(ko[[l]]))
 +
                  names(newind) <- l
 +
                  depind <- c(depind, newind)
 +
                }
 +
              }
 +
            }
 +
          }
 +
        }
 +
      } else {
 +
        dep <- "No dependencies"
 +
      }
 +
      curcols <- colnames(d@output)
 +
      if(!sources) curcols <- curcols[!grepl("Source$", curcols)]
 +
      if(!prevresults) curcols <- curcols[!grepl("Result$", curcols)]
 +
      droploc <- character()
 +
      for(m in curcols) {
 +
        if(!is.numeric(d@output[[m]])) {
 +
          drops <- setdiff(depind[[m]], unique(d@output[[m]]))
 +
          if(length(drops>0)) {
 +
            droploc <- paste(
 +
              droploc,
 +
              paste(
 +
                m,
 +
                paste(drops, collapse = ", "),
 +
                sep = ": "
 +
              ),
 +
              sep = " | "
 +
            )
 +
          }
 +
        }
 +
      }
 +
      if(length(droploc)==0) droploc <- NA
 +
      deptable <- rbind(
 +
        deptable,
 +
        data.frame(
 +
          Ovariable = i,
 +
          Size = nrow(d@output),
 +
          Dependencies = dep,
 +
          Current = paste(curcols, collapse = ", "),
 +
          Dropped = paste(setdiff(names(depind), curcols), collapse = ", "),
 +
          New = paste(setdiff(curcols, names(depind)), collapse = ", "),
 +
          Dropped_locations = droploc
 +
        )
 +
      )
 +
    }
 +
  }
 +
  return(deptable)
 +
}
 +
 
 +
ovashapetest <- function(ova) {
 +
  allr <- rownames(ova@output)
 +
  uniqr <- rownames(unique(ova@output[ova@marginal]))
 +
  cube <- sapply(ova@output[ova@marginal], FUN = function(x) length(unique(x)))
 +
  if(length(allr) == length(uniqr)) {
 +
    cat("All rows have unique marginals.\n")
 +
  } else {
 +
    cat("Warning. All rows do not have unique marginals. Make sure that this is what you want.\n")
 +
  }
 +
  cat("Number of all rows:", length(allr), "\n")
 +
  cat("Number of all rows without Iter: Iter==1", length(ova$Iter[ova$Iter=="1"]),
 +
    "nrow/N", length(allr)/openv$N, "\n")
 +
  cat("Number of unique rows:", length(uniqr), "\n")
 +
  cat("Number of rows in a full array:", prod(cube), "\n")
 +
  oprint(cube)
 +
  nonuniqr <- setdiff(allr, uniqr)
 +
#  cat("Non-unique rows:", nonuniqr, "\n")
 +
#  oprint(head(ova@output[rownames(ova@output) %in% nonuniqr , ]))
 +
  cubesm <- cube[cube>1 & cube<50]
 +
  cubn <- names(cubesm)
 +
  for(i in 2:(length(cubn))) {
 +
    for(j in 1:(i-1)){
 +
      oprint(c(cubn[i], cubn[j]))
 +
      oprint(table(ova@output[[cubn[i]]], ova@output[[cubn[j]]], useNA="ifany"))
 +
    }
 +
  }
 +
 
 +
  for(i in colnames(ova@output)[ova@marginal]) {
 +
    locs <- ova@output[[i]]
 +
    exper <- prod(cube[names(cube) != i])
 +
    oprint(c(i, exper))
 +
    for(j in unique(ova@output[[i]])) {
 +
      cat(j, length(locs[locs == j]), ",") 
 +
    }
 +
  }
 +
}
 +
 
 +
#####################################
 +
# This function can be used to quickly locate indices that do not match between
 +
# two ovariables and thus result in an output with 0 rows.
 +
binoptest <- function(x, y) {
 +
  if(nrow(x@output) == 0) cat(paste("Ovariable", x@name,"has 0 rows in output.\n"))
 +
  if(nrow(y@output) == 0) cat(paste("Ovariable", y@name,"has 0 rows in output.\n"))
 +
  commons <- intersect(colnames(x@output), colnames(y@output))
 +
  commons <- commons[!grepl("Result$", commons)]
 +
  cat("Ovariables have these common columns:\n")
 +
  xt <- x@output
 +
  yt <- y@output
 +
  for (i in commons) {
 +
    cat(i, "with shared locations\n")
 +
    locs <- intersect(x@output[[i[1]]], y@output[[i[1]]])
 +
    if(length(locs)>50) cat(">50 of them\n") else cat(locs, "\n")
 +
    xt <- xt[xt[[i]] %in% locs , ]
 +
    yt <- yt[yt[[i]] %in% locs , ]
 +
    cat("Rows remaining", x@name, nrow(xt), y@name, nrow(yt), "\n")
 +
  }
 +
}
 +
 
 +
#### showLoctable lists locations of each index in the evaluated ovariables in the global environment.
 +
 
 +
showLoctable <- function(name = ".GlobalEnv") {
 +
  loctable <- data.frame()
 +
 
 +
  for(i in ls(name = name)) {
 +
    if(class(get(i)) == "ovariable") {
 +
      for(j in colnames(get(i)@output)) {
 +
        if(!(grepl("Source", j) | grepl("Result", j))) {
 +
          loctable <- rbind(
 +
            loctable,
 +
            data.frame(
 +
              Ovariable = i,
 +
              Index = j,
 +
              Class = paste(class(get(i)@output[[j]]), collapse=" "),
 +
              Marginal = j %in% colnames(get(i)@output)[get(i)@marginal],
 +
              NumLoc = length(unique(get(i)@output[[j]])),
 +
              Locations = paste(head(unique(get(i)@output[[j]])), collapse = " ")
 +
            )
 +
          )
 +
        }
 +
      }
 +
    }
 +
  }
 +
  return(loctable)
 +
}
 +
 
 +
objects.store(showind, binoptest, showLoctable, ovashapetest)
 +
cat("Functions showind, binoptest, showLoctable, ovashapetest stored.\n")
 +
</rcode>
 +
 
 +
'''Functions for Webropol data
 +
 
 +
<rcode name="webropol" embed=1>
 +
#This is code Op_en6007/webropol on page [OpasnetUtils/Drafts]]
 +
 
 +
library(OpasnetUtils)
 +
 
 +
### webropol.convert converts a csv file from Webropol into a useful data.frame.
 +
 
 +
webropol.convert <- function(
 +
  data, # Data.frame created from a Webropol csv file. The first row should contain headings.
 +
  rowfact, # Row number where the factor levels start (in practice, last row + 3)
 +
  textmark = "Other open" # The text that is shown in the heading if there is an open sub-question.
 +
) {
 +
  out <- dropall(data[2:(rowfact - 3) , ])
 +
  subquestion <- t(data[1 , ])
 +
  subquestion <- gsub("\xa0", " ", subquestion)
 +
  subquestion <- gsub("\xb4", " ", subquestion)
 +
  subquestion <- gsub("\n", " ", subquestion)
 +
  #  subquestion <- gsub("\\(", " ", subquestion)
 +
  #  subquestion <- gsub("\\)", " ", subquestion)
 +
  textfield <- regexpr(textmark, subquestion) != -1
 +
  subquestion <- strsplit(subquestion, ":") # Divide the heading into a main question and a subquestion.
 +
  subqtest <- 0 # The previous question name.
 +
  for(i in 1:ncol(out)) {
 +
    #print(i)
 +
    if(subquestion[[i]][1] != subqtest) { # If part of previous question, use previous fact.
 +
      fact <- as.character(data[rowfact:nrow(data) , i]) # Create factor levels from the end of Webropol file.
 +
      fact <- fact[fact != ""] # Remove empty rows
 +
      fact <- gsub("\xa0", " ", fact)
 +
      fact <- gsub("\xb4", " ", fact)
 +
      fact <- gsub("\n", " ", fact)
 +
      fact <- strsplit(fact, " = ") # Separate value (level) and interpretation (label)
 +
    }
 +
    if(length(fact) != 0 & !textfield[i]) { # Do this only if the column is not a text type column.
 +
      out[[i]] <- factor(
 +
        out[[i]],
 +
        levels = unlist(lapply(fact, function(x) x[1])),
 +
        labels = unlist(lapply(fact, function(x) x[2])),
 +
        ordered = TRUE
 +
      )
 +
    }
 +
    subqtest <- subquestion[[i]][1]
 +
  }
 +
  return(out)
 +
}
 +
 
 +
# merge.questions takes a multiple checkbox question and merges that into a single factor.
 +
# First levels in levs have priority over others, if several levels apply to a row.
 +
 
 +
merge.questions <- function(
 +
  dat, # data.frame with questionnaire data
 +
  cols, # list of vectors of column names or numbers to be merged into one level in the factor
 +
  levs, # vector (with the same length as cols) of levels of factors into which questions are merged.
 +
  name # text string for the name of the new factor column in the data.
 +
) {
 +
  for(i in length(cols):1) {
 +
    temp <- FALSE
 +
    for(j in rev(cols[[i]])) {
 +
      temp <- temp | !is.na(dat[[j]])
 +
    }
 +
    dat[[name]][temp] <- levs[i]
 +
  }
 +
  dat[[name]] <- factor(dat[[name]], levels = levs, ordered = TRUE)
 +
  return(dat)
 +
}
 +
 
 +
 
 +
objects.store(webropol.convert, merge.questions)
 +
cat("Functions webropol.convert, merge.questions stored.\n")
 +
</rcode>
 +
 
 +
'''Miscellaneous functions
 +
 
 +
<rcode name="miscellaneous" embed=1>
 +
#This is code Op_en6007/webropol on page [OpasnetUtils/Drafts]]
  
 
library(OpasnetUtils)
 
library(OpasnetUtils)
Line 58: Line 312:
 
}
 
}
 
return(cat("Decisions were forgotten.\n"))
 
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 = NULL, # Base size for graph font (object BS or 18 is used if not given)
 
−
turnx = FALSE, # Turn x axis labels vertically?
 
−
binwidth = NULL # Width of bins
 
−
)
 
−
{
 
−
if(is.null(base_size)) if(exists("BS")) base_size <- BS else base_size <- 18
 
−
if(exists("translate")) {
 
−
ova <- translate(ova)
 
−
x <- translate(x)
 
−
weight <- translate(weight)
 
−
fill <- translate(fill)
 
−
}
 
−
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)
 
 
}
 
}
  
Line 126: Line 335:
 
#DateTime <- as.POSIXct(paste(temperature$Date, temperature$Time), format="%Y-%m-%d %H:%M:%S")
 
#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")
 
#DateHours <- as.POSIXct(paste(price$Date, price$Hours), format="%Y-%m-%d %H:%M:%S")
−
 
−
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 takes a data.frame and fills the cells with NA with each level in that column.
Line 205: Line 378:
 
 
 
return(variable)
 
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 or a data.frame into a SpatialPointsDataFrame.
 
−
  dat, # An evaluated ovariable or data.frame that has coordinate indices.
 
−
  coords = c("LO", "LA"), # The names of the coordinate indices as a character vector, first x then y.
 
−
  proj4string = NULL # Projection identifier or specification as character string. See http://spatialreference.org/
 
−
  # If proj4string is NULL, longitude-latitude system is assumed.
 
−
) {
 
−
  if(class(dat) == "ovariable") temp <- dat@output else
 
−
    if(is.data.frame(dat)) temp <- dat else
 
−
      stop("object must be either evaluated ovariable or data.frame\n")
 
−
 
 
−
  # Transform coordinates into numeric format.
 
−
 
 
−
  for(i in coords) {
 
−
    temp[[i]] <- as.numeric(as.character(temp[[i]]))
 
−
  }
 
−
 
 
−
  # Define the coordinate points first, then add other ovariable output to it.
 
−
 
 
−
  if(is.null(proj4string)) {
 
−
    sp <- SpatialPoints(temp[coords], CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
 
−
    } else {
 
−
      sp <- SpatialPoints(temp[coords], CRS(proj4string))
 
−
    }
 
−
  out <- SpatialPointsDataFrame(sp, temp[!colnames(temp) %in% coords])
 
−
 
 
−
  #Transform the projection to longitude-latitude system.
 
−
  if(!is.null(proj4string)) {
 
−
    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;")
 
 
}
 
}
  
Line 849: Line 747:
 
}
 
}
  
−
orbind2 <- function( # Like orbind but the value is an ovariable.
+
rm(wiki_username)
−
o1, # ovariable whose slots are used in the value.
+
objects.store(list = ls())
−
o2, # ovariable
+
cat("All objects in the global namespace were stored:", ls(), "\n")
−
use_fillna = FALSE, # Do we use fillna to fill in the NA values in indices?
+
 
−
warn = "" # What warning is given if fillna is used?
+
</rcode>
 +
 
 +
'''Functions for GIS data
 +
 
 +
<rcode name="gis" embed=1>
 +
#This is code Op_en6007/gis on page [OpasnetUtils/Drafts]]
 +
 
 +
library(OpasnetUtils)
 +
 
 +
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.
 
) {
 
) {
−
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) {
+
cat("This function MyPointKML is depreciated. Use google.point_kml in OpasnetUtilsExt instead.\n")
−
b <- character()
+
 
−
for(i in colnames(out@output)[out@marginal]) {if(any(is.na(out@output[[i]]))) b <- c(b, i)}
+
    if (is.null(obj))
−
if(length(b) > 0) {
+
        return(list(header = c("<?xml version=\"1.0\" encoding=\"UTF-8\"?>",
−
out@output <- fillna(out@output, b)
+
            "<kml xmlns=\"http://earth.google.com/kml/2.2\">",
−
warning(warn, "\nMissing values had to be filled by function fillna in indices: ", b, "\n")
+
            "<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 or a data.frame into a SpatialPointsDataFrame.
 +
  dat, # An evaluated ovariable or data.frame that has coordinate indices.
 +
  coords = c("LO", "LA"), # The names of the coordinate indices as a character vector, first x then y.
 +
  proj4string = NULL # Projection identifier or specification as character string. See http://spatialreference.org/
 +
  # If proj4string is NULL, longitude-latitude system is assumed.
 +
) {
 +
  if(class(dat) == "ovariable") temp <- dat@output else
 +
    if(is.data.frame(dat)) temp <- dat else
 +
      stop("object must be either evaluated ovariable or data.frame\n")
 +
 
 +
  # Transform coordinates into numeric format.
 +
 
 +
  for(i in coords) {
 +
    temp[[i]] <- as.numeric(as.character(temp[[i]]))
 +
  }
 +
 
 +
  # Define the coordinate points first, then add other ovariable output to it.
 +
 
 +
  if(is.null(proj4string)) {
 +
    sp <- SpatialPoints(temp[coords], CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
 +
    } else {
 +
      sp <- SpatialPoints(temp[coords], CRS(proj4string))
 +
    }
 +
  out <- SpatialPointsDataFrame(sp, temp[!colnames(temp) %in% coords])
 +
 
 +
  #Transform the projection to longitude-latitude system.
 +
  if(!is.null(proj4string)) {
 +
    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")
 
}
 
}
  
−
colnames(out@output)[colnames(out@output) == "Result"] <- paste(o1@name, "Result", sep = "")
+
MyMap <- GetMap(
−
out@marginal <- colnames(out@output) %in% c(xmarg, ymarg)
+
center = c(lat.center, lon.center),
 +
zoom = zoom,
 +
size = size,
 +
RETURNIMAGE = RETURNIMAGE,  
 +
GRAYSCALE = GRAYSCALE,  
 +
verbose = verbose,
 +
...
 +
)
 +
 
 +
PlotOnStaticMap(MyMap) # Plot an empty map.
  
−
return(out)
+
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"
 +
)
 
}
 
}
  
−
rm(wiki_username)
+
MyPlotKML <- function(
−
objects.store(list = ls())
+
shp, # a SpatialPointDataFrame object.
−
cat("All objects in the global namespace were stored:", ls(), "\n")
+
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;")
 +
}
  
 +
objects.store(MyPointKML, ova2spat, MyRmap, MyPlotKML)
 +
cat("Functions MyPointKML, ova2spat, MyRmap, MyPlotKML stored.\n")
 
</rcode>
 
</rcode>
  

Latest revision as of 19:33, 6 June 2017



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") # Old version that fetches all objects, depreciated and not updated.
objects.latest("Op_en6007", code_name = "diagnostics") # Functions for ovariable and model diagnostics: ovashapetest, showLoctable, binoptest
objects.latest("Op_en6007", code_name = "webropol") # Functions for operating with Webropol data
objects.latest("Op_en6007", code_name = "miscellaneous") # Functions for various tasks
objects.latest("Op_en6007", code_name = "gis") # Functions for ovariable, KML and Googl maps interactions

Rationale

Calculations

Functions for ovariable diagnostics showind has problems with get() but this version of code was acceptable [1].

+ Show code

Functions for Webropol data

+ Show code

Miscellaneous functions

+ Show code

Functions for GIS data

+ Show code

See also

References


Related files

<mfanonymousfilelist></mfanonymousfilelist>

Retrieved from "http://www.opasnet.org/testiwiki/index.php?title=OpasnetUtils/Drafts&oldid=40898"