Opasnet map: Difference between revisions

Latest revision as of 11:29, 1 November 2016

This page is a knowledge crystal of subtype method. 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

+ 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

Latest revision as of 11:29, 1 November 2016

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

@@ Line 1: / Line 1: @@
 [[Category:Opasnet]]
 [[Category:Code under inspection]]
-{{variable|moderator=Jouni|stub=Yes}}
+[[Category:R tool]]
+[[Category:Modelling]]
+{{method|moderator=Jouni|stub=Yes}}
 == Question ==
@@ Line 24: / Line 26: @@
 Key guidance:
+* [[OpasnetUtils/Drafts]] contains many functions needed.
+* [[Help:Drawing graphs]] contains guidance for making non-map graphs.
 * [http://cran.r-project.org/web/packages/rgdal/rgdal.pdf rgdal package for R]
 ** [http://www.gdal.org/ GDAL - Geospatial Data Abstraction Library]
@@ Line 32: / Line 36: @@
 == Rationale ==
-All pieces of the puzzle exist already.
+===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
+}}
-=== Dependencies ===
+{{#display_map:
+.928986, 28.706245;
+.937734, 28.712769
+}}
+[http://semantic-mediawiki.org/wiki/Special:MapEditor Map editor and instructions]
 === Calculations ===
 ==== Kuopio buildings with Google pin map ====
+{{attack|#|In function google.point_kml, parameter name is not used for anything!|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 08:42, 26 December 2013 (EET)}}
+{{defend|#|What we need is a function that eats SpatialPointDataFrames and gives KML files with the data points in different formats (pins, coloured marks etc.|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 08:42, 26 December 2013 (EET)}}
 <rcode name='kuorakonmaps'>
@@ Line 50: / Line 70: @@
 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') # Read building data
-shp<-readOGR('PG:host=localhost user=postgres dbname=spatial_db','kuopio_house')
 shp <- shp[1:100 , ] # Save only 100 buildings for demonstration
@@ Line 58: / Line 76: @@
 plotvar	<- shp@data$ika
 nclr	<- 8
-plotclr	<- brewer.pal(nclr,"BuPu")
+plotclr	<- brewer.pal(nclr, "BuPu")
-class	<- classIntervals(plotvar,nclr,style="quantile")
+class	<- classIntervals(plotvar, nclr, style = "quantile")
-colcode	<- findColours(class,plotclr)
 epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
@@ Line 66: / Line 83: @@
 shp2	<- spTransform(shp, epsg4326String)
-kmlname	<- "Kuopio house data"
+dat <- google.point_kml(
-kmldescription <- "Random stuff about here"
+	shp2,
-icon	<- "http://maps.google.com/mapfiles/kml/pal2/icon18.png"
+	kmlname = "Kuopio house data",
-name	<- paste("ika value: ", shp2$ika)
+	kmldescription = "Random stuff about here",
-description <- paste("<b>Value:</b>",shp2$ika,"<br><b>Description:</b>",shp2$kayttotark)
+	name = paste("ika value: ", shp2$ika),
+	description = paste("Age:", shp2$ika, "<br>Description:", shp2$kayttotark),
-dat <- MyPointKML(shp2, kmlname, kmldescription, name, description, icon, colcode)
+	icon = "http://maps.google.com/mapfiles/kml/pal2/icon18.png",
+	col = findColours(class,plotclr)
+)
 google.show_kml_data_on_maps(dat)
@@ Line 78: / Line 97: @@
 </rcode>
+==== GoogleMaps Sorvi MML ====
-==== GoogleMaps Sorvi MML TEST ====
-{{attack|#|Something wrong with sorvi.|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 12:54, 23 December 2013 (EET)}}
 <rcode name='gmapspsqltest3'>
@@ Line 89: / Line 105: @@
 library(sorvi)
 library(rgdal)
+library(maptools)
-data(MML)
+shp <- LoadData("kuntarajat.maa.shp")
-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")
@@ Line 99: / Line 114: @@
 shp2<-spTransform(shp,epsg4326String)
-out<-sapply(slot(shp2,"polygons"),function(x
+out<-sapply(slot(shp2,"polygons"),function(x)
 	{kmlPolygon(
 		x,
@@ Line 210: / Line 225: @@
 </rcode>
-==== Static GoogleMaps test ====
+====Rasters on Google maps====
+<rcode graphics=1>
+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(
+.942533,
+.945662,
+.931893,
+.925618,
+.936103,
+.945605,
+.986605,
+.992419,
+.004996,
+.044465,
+.950091,
+.992843,
+.968536,
+.012646,
+.018907,
+.935284,
+.003882
+),
+LO = c(
+.936602,
+.971964,
+.907677,
+.929521,
+.926474,
+.953424,
+.971706,
+.983251,
+.002756,
+.104916,
+.984538,
+.984581,
+.953897,
+.001447,
+.022261,
+.937717,
+.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"
+	)
+</rcode>
+==== Static GoogleMaps ====
 <rcode name='staticmapstest' graphics=1 embed=1>
@@ Line 239: / Line 422: @@
 * [https://gist.github.com/2292662 Plotting sp-objects to Google Maps with different CRSs]
+* [http://semantic-mediawiki.org/wiki/Special:MapEditor Map editor]
+* [http://www.findlatitudeandlongitude.com/batch-geocode/ Batch geocoding]: give addresses, the website will provide geolocations. Works for a few hundred addresses.
 ==Keywords==