Opasnet map: Difference between revisions

Revision as of 09:31, 29 January 2014

This page is a knowledge crystal of subtype variable. The page identifier is Op_en5490
Moderator:Jouni (see all)
This page is a stub. You may improve it into a full page.
Upload data {{#opasnet_base_link:Op_en5490}}

Question

How should GIS data be handled and visualised in Opasnet?

Answer

Original GIS data is stored in a PostgreSQL database.
Data is accessed by R using the RGDAL package and ?? connection.
Data in manipulated in R.
Ovariables are converted to SpatialPointsDataFrame objects using ova2spat function.
Data is displayed by producing a KML file with ?? package.
- The KML file is created with MyPointKML function, if pins are shown.
- The KML file is saved at the R-tools server.
- Google Maps is used to show the KML file on a web page.

Currently, ova2spat and MyPointKML functions are located in OpasnetUtils/Drafts. Therefore, you need this command:

objects.latest("Op_en6007", code_name = "answer") # OpasnetUtils/Drafts. We need MyPointsKML and ova2spat.

This is the projection that the National Land Survey Finland uses: [1].

Key guidance:

OpasnetUtils/Drafts contains many functions needed.
Help:Drawing graphs contains guidance for making non-map graphs.
rgdal package for R
- GDAL - Geospatial Data Abstraction Library
sp package for R: contains descriptions about the key spatial objects.
Changing the coordinate system of Spatial objects

Rationale

Wiki maps

{{#display_map: 62.91826544557644,28.69525909423828~ ~Valuma-alue; 62.92647012476198,28.637065887451172 |rectangles=62.94068612542196,28.726844787597656:62.92373548704899,28.702983856201172~Luikonlahden kaivos~ ~ ~ ~ ~#F2E1E1 |zoom=10 |centre=62.91826544557644,28.69525909423828 }}

Map editor and instructions

Calculations

Kuopio buildings with Google pin map

⇤--#: . In function google.point_kml, parameter name is not used for anything! --Jouni (talk) 08:42, 26 December 2013 (EET) (type: truth; paradigms: science: attack)

←--#: . What we need is a function that eats SpatialPointDataFrames and gives KML files with the data points in different formats (pins, coloured marks etc. --Jouni (talk) 08:42, 26 December 2013 (EET) (type: truth; paradigms: science: defence)

+ Show code - Hide code


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

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

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

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

dat <- google.point_kml(
	shp2, 
	kmlname = "Kuopio house data", 
	kmldescription = "Random stuff about here", 
	name = paste("ika value: ", shp2$ika), 
	description = paste("Age:", shp2$ika, "<br>Description:", shp2$kayttotark), 
	icon = "http://maps.google.com/mapfiles/kml/pal2/icon18.png", 
	col = findColours(class,plotclr)
)

google.show_kml_data_on_maps(dat)

GoogleMaps Sorvi MML

⇤--#: . Something wrong with sorvi. --Jouni (talk) 12:54, 23 December 2013 (EET) (type: truth; paradigms: science: attack)

+ Show code - Hide code


library(OpasnetUtils)
library(OpasnetUtilsExt)
library(sorvi)
library(rgdal)
library(maptools)

shp <- LoadData("kuntarajat.maa.shp")

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

Google with shapefiles

+ Show code - Hide code


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

Google show data from url on map

+ Show code - Hide code

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

Google circles

+ Show code - Hide code

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)

Rasters on Google maps

+ Show code - Hide code

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

Static GoogleMaps

+ Show code - Hide code


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.

Keywords

References

Related files

Opasnet map: Difference between revisions

Revision as of 09:31, 29 January 2014

Contents

Question

Answer

Rationale

Wiki maps

Calculations

Kuopio buildings with Google pin map

GoogleMaps Sorvi MML

Google with shapefiles

Google show data from url on map

Google circles

Rasters on Google maps

Static GoogleMaps

See also

Keywords

References

Related files

Navigation menu