Difference between revisions of "Opasnet map"

From Testiwiki
Jump to: navigation, search
(first draft)
 
(→‎See also)
 
(20 intermediate revisions by 3 users not shown)
Line 1: Line 1:
 
[[Category:Opasnet]]
 
[[Category:Opasnet]]
−
{{variable|moderator=Jouni|stub=Yes}}
+
[[Category:Code under inspection]]
 +
[[Category:R tool]]
 +
[[Category:Modelling]]
 +
{{method|moderator=Jouni|stub=Yes}}
  
 
== Question ==
 
== Question ==
Line 11: Line 14:
 
* Data is accessed by [[R]] using the RGDAL package and ?? connection.
 
* Data is accessed by [[R]] using the RGDAL package and ?? connection.
 
* Data in manipulated in [[R]].
 
* Data in manipulated in [[R]].
 +
* Ovariables are converted to SpatialPointsDataFrame objects using ova2spat function.
 
* Data is displayed by producing a [[:en:Keyhole markup language|KML]] file with ?? package.
 
* Data is displayed by producing a [[:en:Keyhole markup language|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.
 
** The KML file is saved at the R-tools server.
 
** Google Maps is used to show the KML file on a web page.
 
** 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: [http://spatialreference.org/ref/epsg/3067/].
 +
 +
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]
 +
* [http://cran.r-project.org/web/packages/sp/vignettes/intro_sp.pdf sp package for R]: contains descriptions about the key spatial objects.
 +
* [http://rwiki.sciviews.org/doku.php?id=tips:spatial-data:change_crs Changing the coordinate system of Spatial objects]
 +
  
 
== Rationale ==
 
== 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
 +
}}
 +
 
 +
 
 +
 
 +
{{#display_map:
 +
62.928986, 28.706245;
 +
62.937734, 28.712769
 +
}}
 +
 
 +
[http://semantic-mediawiki.org/wiki/Special:MapEditor Map editor and instructions]
 +
 
 +
=== Calculations ===
  
 +
==== Kuopio buildings with Google pin map ====
  
−
=== Dependencies ===
+
{{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)}}
  
−
=== Formula ===
+
{{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'>
 +
 
 +
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)
 +
 
 +
</rcode>
 +
 
 +
==== GoogleMaps Sorvi MML ====
 +
 
 +
<rcode name='gmapspsqltest3'>
 +
 
 +
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)
 +
 
 +
</rcode>
 +
 
 +
==== Google with shapefiles ====
 +
 
 +
<rcode name='gmapspsqltest2'>
 +
 
 +
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()
 +
 
 +
</rcode>
 +
 
 +
==== Google show data from url on map ====
 +
 
 +
<rcode name='gutilstest'>
 +
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')
 +
</rcode>
 +
 
 +
==== Google circles ====
 +
 
 +
<rcode name='gmaptest'>
 +
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)
 +
 
 +
</rcode>
 +
 
 +
====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(
 +
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"
 +
)
 +
 
 +
</rcode>
 +
 
 +
==== Static GoogleMaps ====
 +
 
 +
<rcode name='staticmapstest' graphics=1 embed=1>
 +
 
 +
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.
 +
 
 +
</rcode>
  
 
==See also==
 
==See also==
 +
 +
* [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==
 
==Keywords==

Latest revision as of 11:29, 1 November 2016



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:


Rationale

Wiki maps

Loading map...


Loading map...

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)

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

+ Show code

GoogleMaps Sorvi MML

+ Show code

Google with shapefiles

+ Show code

Google show data from url on map

+ Show code

Google circles

+ Show code

Rasters on Google maps

+ Show code

Static GoogleMaps

+ Show code

See also

Keywords

References


Related files

<mfanonymousfilelist></mfanonymousfilelist>

Retrieved from "http://www.opasnet.org/testiwiki/index.php?title=Opasnet_map&oldid=39994"