Difference between revisions of "Opasnet map"

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

objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]]; we need MyPointsKML and ova2spat.

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

shp <- shp[1:100 , ] # Save only 100 buildings for demonstration

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("ika value: ", shp2$ika)
description <- paste("<b>Value:</b>",shp2$ika,"<br><b>Description:</b>",shp2$kayttotark)

dat <- MyPointKML(shp2, kmlname, kmldescription, name, description, icon, colcode)

google.show_kml_data_on_maps(dat)

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

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=''
)

cat("PostgreSQL test with coloured data areas.\n")

google.show_kml_data_on_maps(data)

cat("PostgreSQL Test with default data\n")

google.show_data_on_maps()

library(OpasnetUtils)
library(OpasnetUtilsExt)
google.show_kml_on_maps(url='http://api.flickr.com/services/feeds/geo/?g=322338@N20&lang=en-us&format=feed-georss')

library(OpasnetUtilsExt)

data = list()
data[[1]] = c(62.8925, 27.678333, 9683, '#ff0000', 'kuopio') 
data[[2]] = c(65.016667, 25.466667, 14398, '#00ff00', 'oulu') 
data[[3]] = c(60.170833, 24.9375, 59623, '#0000ff', 'helsinki') 

google.circles(data)

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

altiste <- "Ni"

pitoisuus <- new("ovariable",
	name = "pitoisuus",
	data = tidy(op_baseGetData("opasnet_base", "Op_fi3280", series_id = 5509), objname= "pitoisuus")
)

# Remove empty rows.
pitoisuus@data <- pitoisuus@data[pitoisuus@data$pitoisuusResult != "" , ]

# Convert values that are below detection limit to ranges 0 - detection limit.
pitoisuus@data$pitoisuusResult <- gsub("<", "0-", pitoisuus@data$pitoisuusResult)

# Convert different time formats used to POSIXct time format.
temp <- as.POSIXct(pitoisuus@data$Aika, format = "%d.%m.%Y %H:%M")
temp <- ifelse(is.na(temp), as.POSIXct(pitoisuus@data$Aika, format = "%d.%m.%Y %H.%M"), temp)
temp <- ifelse(is.na(temp), as.POSIXct(pitoisuus@data$Aika, format = "%d.%m.%Y"), temp)
temp <- as.POSIXct(temp, origin = "1970-01-01 02:00.00 UTC")
pitoisuus@data$Aika <- temp

pitoisuus <- EvalOutput(pitoisuus, N = 1)

data <- pitoisuus@output[pitoisuus@output$Altiste == altiste, ]
data.time <- data

########### Concentrations on map

locs <- data.frame(Paikka = c(
	"Lumijärvi",
	"Ylä-Lumijärvi",
	"Kivijärvi",
	"Lumijoki, Kivijärveä lähin piste",
	"Lumijoki, keskimmäinen piste",
	"Lumijoki, Lumijärveä lähin piste",
	"Lampi kaivosalueen reunalla 1",
	"Lampi kaivosalueen reunalla 2",
	"Puro Kalliojärven ja Salmisen välillä",
	"Tuhkajoki",
	"Kortelampi",
	"Kärsälampi",
	"Rengaskaivo, Pappila",
	"Kalliojärvi",
	"Kalliojoki",
	"Lumijoki",
	"Salmisenpuro"
),
LA = c(
	63.942533,
	63.945662,
	63.931893,
	63.925618,
	63.936103,
	63.945605,
	63.986605,
	63.992419,
	64.004996,
	64.044465,
	63.950091,
	63.992843,
	63.968536,
	64.012646,
	64.018907,
	63.935284,
	64.003882
),
LO = c(
	27.936602,
	27.971964,
	27.907677,
	27.929521,
	27.926474,
	27.953424,
	27.971706,
	27.983251,
	28.002756,
	28.104916,
	27.984538,
	27.984581,
	27.953897,
	28.001447,
	28.022261,
	27.937717,
	27.999601
))

oprint(locs)

data <- merge(data, locs)
result <- "pitoisuusResult"

title <- paste(altiste, "Talvivaaran läheisyydessä (", data$pitoisuusUnit[1], ")", sep = "")
coordinates(data)=c("LO","LA")

# Plot the data
proj4string(data)<-("+init=epsg:4326")
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
shp<-spTransform(data,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) <- 32
nrow(rast) <- 32

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

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

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;'>", title, "</span>\n")

colorstrip(colors, steps)

#Plot data

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

ggplot(data.time, aes(x = Paikka, y = pitoisuusResult, colour = pitoisuusResult)) + 
	geom_point(size = 5) + # point size
	theme_grey(base_size = 24) + # font etc. http://sape.inf.usi.ch/quick-reference/ggplot2/themes
	labs(
		title = paste(altiste, "Talvivaaran läheisyydessä"), 
		y = paste("Pitoisuus,", data$pitoisuusUnit[1]),
		x = "Mittauspaikka"
	) 

########### Concentrations in time

ggplot(data.time, aes(x = Aika, y = pitoisuusResult, colour = Paikka)) + 
	geom_point(size = 5) + 
	theme_grey(base_size = 24) +
	labs(
		title = paste(altiste, "Talvivaaran läheisyydessä"), 
		y = paste("Pitoisuus,", data$pitoisuusUnit[1]),
		x = "Aika"
	) 

library(RgoogleMaps)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(OpasnetUtils)
library(OpasnetUtilsExt)

objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]]. We need MyRmap.

shp	<- readOGR('PG:host=localhost user=postgres dbname=spatial_db','kuopio_house')
proj4string(shp) <- ("+init=epsg:3067") # The map projection of NLS Finland.

#oprint(shp@data[1:100 , ])

epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
shp2	<- spTransform(shp,epsg4326String) # Convert to longitude-latitude projection.

MyRmap(shp2, plotvar = "ika", legend_title = "Ikä", numbins = 8, pch = 19, cex = 0.3) # Draw the map.