Difference between revisions of "Wikisym 2013 Demo"

a <-3
b <-4
a +b

sum <- num1 + num2
mul <- num1 * num2
cat(paste("<span style='font-size: 1.5em;font-weight:bold;'>Sum: ",sum,"</span>\n",sep=''))
cat(paste("<span style='font-size: 1.5em;font-weight:bold;'>Mul: ",mul,"</span>\n",sep=''))

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

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(OpasnetUtils)
library(OpasnetUtilsExt)
library(rgdal)
library(maptools)
library(raster)

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)
}

shp_o <- spatial_db_query(paste('SELECT * FROM metsapalot_pm2_5_europe_2006 WHERE pm2_5 > ',treshold,' AND kk = ',kk,' AND pv = ',pv,';',sep=''))

#THIS IS A KEY PART OF FILTER CODE: As you are calling in a Query you need to tell R what the Lat Lon attributes are of the data
coordinates(shp_o)=c("lon","lat")

#plotvar removed as points all same value -- in future we may need an 'IF' clause here
#all reference to colcode removed in this example

 
#proj <- CRS("+ellps=sphere +a=127.4 +e=0 +proj=stere +lat_0=90 +lon_0=-32 +lat_ts=60 +x_0=8 +y_0=110")
proj <- ("+init=epsg:4326")
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
proj4string(shp_o) <- proj
shp<-spTransform(shp_o,epsg4326String)

#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)<-max(shp$emep_i)*2
nrow(rast)<-max(shp$emep_j)*2

cat("<span style='font-size: 1.2em;font-weight:bold;'>Wild Fire PM 2.5 Emissions in Europe</span>\n")

#Rasterize point data
rast2<-rasterize(shp, rast, shp$pm2_5, fun=mean)

#max(shp$pm2_5)
#min(shp$pm2_5)

#colors <- rainbow(11,start=0,end=0.15)

steps <- round(min(shp$pm2_5)):round(max(shp$pm2_5))
colors <- rainbow(length(steps), start=0, end=0.50)

colorstrip(colors, steps)

#Plot data
google.show_raster_on_maps(rast2, col=colors, style="height:300px; width:300px;")

library("OpasnetUtils")
library("xtable")

# From opasnet
thedata <- opbase.data("Op_en2949")
print(xtable(thedata),type="html")