Difference between revisions of "Sandbox"

library(OpasnetUtils)
a <- new("ovariable", output = data.frame(dummy=NA))
b <- new("ovariable", output = data.frame(a=1:4))

merge(a,b)

#code goes here
library(RgoogleMaps)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(OpasnetBaseUtils)

shp<-readOGR('PG:host=localhost user=postgres dbname=spatial_db','kuopio_house')
plotvar<-shp@data$ika
nclr<-8
plotclr<-brewer.pal(nclr,"Spectral")
class<-classIntervals(plotvar,nclr,style="quantile")
colcode<-findColours(class,plotclr)

epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
proj4string(shp)<-("+init=epsg:3067")
shp2<-spTransform(shp,epsg4326String)
#get marker information first 10 points
mymarkers<-cbind.data.frame(lat=c(shp2@coords[,2]),lon=c(shp2@coords[,1]),color=colcode);

#get the bounding box:
bb <- qbbox(lat = mymarkers[,"lat"], lon = mymarkers[,"lon"])

#MyMap function without the "file destination" parameter
MyRmap<-function (lonR, latR, center, size = c(640, 640),  
    MINIMUMSIZE = FALSE, RETURNIMAGE = TRUE, GRAYSCALE = FALSE, 
    NEWMAP = TRUE, zoom, verbose = 1, ...) 
{
    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)
        ur <- LatLon2XY(latR[2], lonR[2], zoom)
        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")
    }
    return(GetMap(center = c(lat.center, lon.center), zoom = zoom, 
        size = size, RETURNIMAGE = RETURNIMAGE, 
        GRAYSCALE = GRAYSCALE, verbose = verbose, ...))
}

MyMap<-MyRmap(bb$lonR,bb$latR,maptype="mobile")

PlotOnStaticMap(MyMap)

PlotOnStaticMap(MyMap,lat=mymarkers[,"lat"],lon=mymarkers[,"lon"],pch=19,cex=0.3,col=colcode,add=T)

legend("topleft", legend=names(attr(colcode, "table")),title="Ika", fill=attr(colcode, "palette"),  cex=1.0, bty="y",bg="white")

library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(OpasnetBaseUtils)
library(RODBC)

shp <- spatial_db_query(paste('SELECT * FROM kuopio_house WHERE ika >= ',age,';',sep=''))

coordinates(shp)=c("y_koord","x_koord")

#shp<-readOGR('PG:host=localhost user=postgres dbname=spatial_db','kuopio_house')
plotvar<-shp@data$ika
nclr<-8
plotclr<-brewer.pal(nclr,"BuPu")
class<-classIntervals(plotvar,nclr,style="quantile")
colcode<-findColours(class,plotclr)

epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
proj4string(shp)<-("+init=epsg:3067")
shp2<-spTransform(shp,epsg4326String)


kmlname<-"Kuopio house data"
kmldescription<-"Random stuff about here"
icon<-"http://maps.google.com/mapfiles/kml/pal2/icon18.png"
name<-paste("Value: ",shp2$ika)
description <- paste("<b>Ikä:</b>",shp2$ika,"<br><b>Rakennustunnus:</b>",shp2$rakennustunnus)

MyPointKML<-function(obj = NULL, kmlname = "", kmldescription = "", name = NULL, description = "", icon = "http://maps.google.com/mapfiles/kml/pal4/icon24.png",col=NULL) 
{
    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))
    }

    col2kmlcolor <- function(col) paste(rev(sapply(col2rgb(col, TRUE), function(x) sprintf("%02x", x))), collapse = "")
    kml <- kmlStyle <- ""
    kmlHeader <- c("<?xml version=\"1.0\" encoding=\"UTF-8\"?>","<kml xmlns=\"http://earth.google.com/kml/2.2\">", "<Document>")
    kmlFooter <- c("</Document>", "</kml>")
    for (i in 1:nrow(obj)) {
        point <- obj[i, ]
        pt_name = name[i]
        pt_description = description[i]
        pt_style <- paste("#style", ifelse(length(icon) == 1, 1, i), sep = "")
        kml <- append(kml, "<Placemark>")
        kml <- append(kml, paste("  <description><![CDATA[",pt_description, "]]></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 = ""))
    
}



data <- MyPointKML(shp2,kmlname,kmldescription,name,description,icon,colcode)
google.show_kml_data_on_maps(data)

library(OpasnetBaseUtils)
library(sorvi)
library(rgdal)

data(MML)

shp <- MML[["1_milj_Shape_etrs_shape"]][["kunta1_p"]]

#epsg3857String <- CRS("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
proj4string(shp)<-("+init=epsg:3047")
shp2<-spTransform(shp,epsg4326String)

out<-sapply(slot(shp2,"polygons"),function(x){kmlPolygon(x,name="nimi",col='#df0000aa',lwd=1,border='black',description="selite") })

data<-paste(
paste(kmlPolygon(kmlname="This will be layer name", kmldescription="<i>More info about layer here</i>")$header, collapse="\n"),
paste(unlist(out["style",]), collapse="\n"),
paste(unlist(out["content",]), collapse="\n"),
paste(kmlPolygon()$footer, collapse="\n"),
sep=''
)

google.show_kml_data_on_maps(data)

library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(OpasnetBaseUtils)

shp<-readOGR('PG:host=localhost user=postgres dbname=spatial_db','watson_wkt')

plotvar<-shp@data$value_inhalation
nclr<-8
plotclr<-brewer.pal(nclr,"BuPu")
class<-classIntervals(plotvar,nclr,style="quantile")
colcode<-findColours(class,plotclr)
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
proj4string(shp)<-("+init=epsg:3035")
shp2<-spTransform(shp,epsg4326String)


out<-sapply(slot(shp2,"polygons"),function(x){kmlPolygon(x,name=as(shp2,"data.frame")[slot(x,"ID"),"country_code"],col=colcode[[((as.numeric(slot(x,"ID"))+1))]],lwd=1,border='black',description=paste("Value:",as(shp2,"data.frame")[slot(x,"ID"),"value_inhalation"])) })

data<-paste(
paste(kmlPolygon(kmlname="This will be layer name", kmldescription="<i>More info about layer here</i>")$header, collapse="\n"),
paste(unlist(out["style",]), collapse="\n"),
paste(unlist(out["content",]), collapse="\n"),
paste(kmlPolygon()$footer, collapse="\n"),
sep=''
)

google.show_kml_data_on_maps(data)

library('OpasnetBaseUtils')
cat("<span style='font-size: 1.2em;font-weight:bold;'>PostgreSQL Test</span>\n")
google.show_data_on_maps()


google.show_data_on_maps(table='kuopio_house',database='spatial_db',fields=c('ika','ika','the_geom'))

library(OpasnetBaseUtils)

csv <- opasnet.csv("2/25/Russian_elections_2011_results.csv")

print(csv[1:10,1:10])

library(OpasnetBaseUtils)
pumps.model <- opasnet.data('c/cc/Test_bugs_model.txt')

library('rbugs')
data(pumps)
pumps.data <- list(t = pumps$t, x = pumps$x, N = nrow(pumps))
inits <- list(alpha = 1, beta = 1)
parameters <- c("theta", "alpha", "beta")
pumps.sim <- bugs.run(data = pumps.data, list(inits), parameters,pumps.model, n.chains = 1, n.iter = 1000)
## MCMC Analysis
library("coda")
pumps.mcmc <- as.mcmc(pumps.sim$chain1)
summary(pumps.mcmc)
effectiveSize(pumps.mcmc)
## End(Not run)

# This code does not work. Should test the one in test wiki but i don't remember the location.
library(rjags)
N <- 1000
x <- rnorm(N, 0, 5)
 
example.bug <- "model {
	for (i in 1:N) {
		x[i] ~ dnorm(mu, tau)
	}
	mu ~ dnorm(0, .0001)
	tau <- pow(sigma, -2)
	sigma ~ dunif(0, 100)
}"

jags.model

jags <- jags.model(example1.bug,
                   data = list('x' = x,
                               'N' = N),
                   n.chains = 4,
                   n.adapt = 100)
 
update(jags, 1000)
 
jags.samples(jags,
             c('mu', 'tau'),
             1000)

cat("Above should be included code\n")

library(OpasnetBaseUtils)
library(ggplot2)
library(xtable)

saanto.siemenet <- op_baseGetData("opasnet_base", "Op_fi2633")[,-c(1,2,7)] # Jatropan siementen saanto viljelystä
saanto.öljy     <- op_baseGetData("opasnet_base", "Op_fi2634")[,-c(1,2,5)] # Öljyn saanto jatropan siemenistä
saanto.diesel   <- op_baseGetData("opasnet_base", "Op_fi2632")[,-c(1,2,5)] # Biodieselin saanto jatropaöljystä
viljelyala      <- op_baseGetData("opasnet_base", "Op_fi2642")[,-c(1,2)] # Jatropan viljelyalueet
päästö.ilmasto  <- op_baseGetData("opasnet_base", "Op_fi2547")[,-c(1,2)] # Jatropan viljelyn ilmastovaikutukset
päästö.sosiaali <- op_baseGetData("opasnet_base", "Op_fi2552")[,-c(1,2)] # Jatropan viljelyn sosiaaliset vaikutukset
päästö.ekosyst  <- op_baseGetData("opasnet_base", "Op_fi2548")[,-c(1,2)] # Jatropan viljelyn ekosysteemivaikutukset
P               <- op_baseGetData("opasnet_base", "Op_fi2539")[,-c(1,2,7)] # Jatropan käyttö bioenergian lähteenä

colnames(saanto.siemenet)[4] <- "siemenet"
colnames(saanto.öljy)[2] <- "öljy"
colnames(saanto.diesel)[2] <- "diesel"
saanto <- merge(saanto.siemenet, saanto.öljy)
saanto <- merge(saanto, saanto.diesel)
saanto[,9] <- saanto$siemenet * saanto$öljy * saanto$diesel * ala
colnames(saanto)[9] <- "saanto (kg/a)"

P <- PTable(P, n)
saanto <- merge(P, saanto)

if(length(divisions)>1) divisions <- as.list(saanto[, divisions]) else divisions <- saanto[, divisions]
out1 <- as.data.frame(as.table(tapply(saanto[, 10], divisions, mean))) 
out1 <- dropall(out1[!is.na(out1$Freq), ])

print(xtable(out1), type = 'html')

out2 <- as.data.frame(as.table(tapply(saanto[, 10], list(saanto[, divisions2], saanto$obs), mean))) 
out2 <- dropall(out2[!is.na(out2$Freq), ])
out2[1:10, ]
ggplot(out2, aes(x = Freq, weight = 1, fill = Var1)) +geom_density() 
## Jostain syystä vain osa kuvista piirtyy oikein, riippuen mitä parametreja valitaan. En ymmärrä syytä.

######################################
## dropall pudottaa data.framesta pois kaikki faktorien sellaiset levelit, joita ei käytetä.
## parametrit: x = data.frame 

dropall <- function(x){
    isFac <- NULL
    for (i in 1:dim(x)[2]){isFac[i] = is.factor(x[ , i])}

    for (i in 1:length(isFac)){
        x[, i] <- x[, i][ , drop = TRUE]
        }
    return(x)
    }
########################################

#########################################
## PTable muuntaa arvioinnin todennäköisyystaulun sopivaan muotoon arviointia varten.
## Parametrit: P = todennäköisyystaulu Opasnet-kannasta kaivettuna.
##             n = iteraatioiden lukumäärä Monte Carlossa
## Todennäköisyystaulun sarakkeiden on oltava: Muuttuja, Selite, Lokaatio, P
## Tuotteena on Monte Carloa varten tehty taulu, jonka sarakkeina ovat
## n (iteraatio) ja kaikki todennäköisyystaulussa olleet selitteet, joiden riveille on arvottu
## lokaatiot niiden todennäköisyyksien mukaisesti, jotka todennäköisyystaulussa oli annettu.

PTable <- function(P, n) {
Pt <- unique(P[,c("Muuttuja", "Selite")])
Pt <- data.frame(Muuttuja = rep(Pt$Muuttuja, n), Selite = rep(Pt$Selite, n), obs = rep(1:n, each = nrow(Pt)), P = runif(n*nrow(Pt), 0, 1))
for(i in 2:nrow(P)){P$Result[i] <- P$Result[i] + ifelse(P$Muuttuja[i] == P$Muuttuja[i-1] & P$Selite[i] == P$Selite[i-1], P$Result[i-1], 0)}
P <- merge(P, Pt)
P <- P[P$P <= P$Result, ]
Pt <- as.data.frame(as.table(tapply(P$Result, as.list(P[, c("Muuttuja", "Selite", "obs")]), min)))
colnames(Pt) <- c("Muuttuja", "Selite", "obs", "Result")
Pt <- Pt[!is.na(P$Result), ]
P <- merge(P, Pt)
P <- P[, !colnames(P) %in% c("Result", "P", "Muuttuja")]
P <- reshape(P, idvar = "obs", timevar = "Selite", v.names = "Lokaatio", direction = "wide")
colnames(P) <- ifelse(substr(colnames(P), 1, 9) == "Lokaatio.", substr(colnames(P), 10,30), colnames(P))
return(P)
}
######################################

# The area of the rect is
width *  height

library(OpasnetUtils)
library(OpasnetUtilsExt)
library(ggplot2)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(raster)

data <- MassHEIS.data()

# Plot the data
coordinates(data)=c("longitude","latitude")
proj4string(data)<-("+init=epsg:4326")
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
shp<-spTransform(data,epsg4326String)


start <- min(shp$anmean)
end <- max(shp$anmean)

rasters = list()
years = c(2000,2001,2002)

for(y in 1:length(years))
{
 #Create blank raster
 rast<-raster()

 s <- shp[(shp$Year == years[y]),]

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

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

 #Rasterize point data
 rasters[[y]] <- rasterize(s, rast, s$anmean, fun=mean)
}

steps <- approx(c(start,end),n=6)$y
colors <- rev(rainbow(length(steps), start=0, end=0.50))

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

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

cat("<span style='font-size: 1.2em;font-weight:bold;'>Massachusetts annual mean PM 2.5 microns</span>\n")

colorstrip(colors, steps)

#Plot data

print(rasters)

s <- stack(rasters)

#names(s) <- c('2010','2009') 

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