Difference between revisions of "Opasnet map"

From Testiwiki
Jump to: navigation, search
(→‎Static GoogleMaps test: function reorganised)
(→‎See also)
 
(12 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 
[[Category:Opasnet]]
 
[[Category:Opasnet]]
 
[[Category:Code under inspection]]
 
[[Category:Code under inspection]]
−
{{variable|moderator=Jouni|stub=Yes}}
+
[[Category:R tool]]
 +
[[Category:Modelling]]
 +
{{method|moderator=Jouni|stub=Yes}}
  
 
== Question ==
 
== Question ==
Line 24: Line 26:
  
 
Key guidance:
 
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://cran.r-project.org/web/packages/rgdal/rgdal.pdf rgdal package for R]
 
** [http://www.gdal.org/ GDAL - Geospatial Data Abstraction Library]
 
** [http://www.gdal.org/ GDAL - Geospatial Data Abstraction Library]
Line 32: Line 36:
 
== 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
 +
}}
  
−
=== Dependencies ===
+
 
 +
 
 +
{{#display_map:
 +
62.928986, 28.706245;
 +
62.937734, 28.712769
 +
}}
 +
 
 +
[http://semantic-mediawiki.org/wiki/Special:MapEditor Map editor and instructions]
  
 
=== Calculations ===
 
=== Calculations ===
  
 
==== Kuopio buildings with Google pin map ====
 
==== 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'>
 
<rcode name='kuorakonmaps'>
Line 50: Line 70:
 
library(OpasnetUtilsExt)
 
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
 
shp <- shp[1:100 , ] # Save only 100 buildings for demonstration
Line 58: Line 76:
 
plotvar <- shp@data$ika
 
plotvar <- shp@data$ika
 
nclr <- 8
 
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")
 
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
Line 66: Line 83:
 
shp2 <- spTransform(shp, epsg4326String)
 
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)
 
google.show_kml_data_on_maps(dat)
Line 78: Line 97:
 
</rcode>
 
</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'>
 
<rcode name='gmapspsqltest3'>
Line 89: Line 105:
 
library(sorvi)
 
library(sorvi)
 
library(rgdal)
 
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")
 
#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)
 
shp2<-spTransform(shp,epsg4326String)
  
−
out<-sapply(slot(shp2,"polygons"),function(x
+
out<-sapply(slot(shp2,"polygons"),function(x)
 
{kmlPolygon(
 
{kmlPolygon(
 
x,
 
x,
Line 210: Line 225:
 
</rcode>
 
</rcode>
  
−
==== Static GoogleMaps test ====
+
====Rasters on Google maps====
  
−
<rcode name='staticmapstest' graphics='1'>
+
<rcode graphics=1>
−
 
+
library(OpasnetUtils)
−
library(RgoogleMaps)
+
library(OpasnetUtilsExt)
 
library(rgdal)
 
library(rgdal)
 
library(maptools)
 
library(maptools)
 
library(RColorBrewer)
 
library(RColorBrewer)
 
library(classInt)
 
library(classInt)
−
library(OpasnetUtils)
+
library(raster)
−
library(OpasnetUtilsExt)
+
 
 +
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)
  
−
shp <- readOGR('PG:host=localhost user=postgres dbname=spatial_db','kuopio_house')
+
data <- merge(data, locs)
−
#oprint(shp@data[1:100 , ])
+
result <- "pitoisuusResult"
  
−
plotvar <- shp@data$ika
+
title <- paste(altiste, "Talvivaaran läheisyydessä (", data$pitoisuusUnit[1], ")", sep = "")
−
nclr <- 8
+
coordinates(data)=c("LO","LA")
  
 +
# Plot the data
 +
proj4string(data)<-("+init=epsg:4326")
 
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
 
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
−
proj4string(shp) <- ("+init=epsg:3067")
+
shp<-spTransform(data,epsg4326String)
−
shp2 <- spTransform(shp,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)
  
−
# MyMap function without the "file destination" parameter
+
#Plot data
−
# Requires RgoogleMaps package
 
  
−
MyRmap <- function (
+
google.show_raster_on_maps(rast2, col=colors, style="height:500px;")
−
shp2, # a spatial data object
 
−
# lonR, # longitude of data points
 
−
# latR, # latitude of data points
 
−
center, # center of the map
 
−
size = c(640, 480), # size of the map. This produces the right dimensions in Opasnet.
 
−
MINIMUMSIZE = FALSE,
 
−
RETURNIMAGE = TRUE,
 
−
GRAYSCALE = FALSE,
 
−
NEWMAP = TRUE,
 
−
zoom,
 
−
verbose = 1,
 
−
pch = 19,
 
−
cex = 0.3,
 
−
legend_title = "",
 
−
...
 
−
) {
 
−
plotclr <- brewer.pal(nclr,"Spectral")
 
−
class <- classIntervals(plotvar,nclr,style="quantile")
 
−
colcode <- findColours(class,plotclr)
 
  
−
latR <- shp2@coords[ , 2]
+
ggplot(data.time, aes(x = Paikka, y = pitoisuusResult, colour = pitoisuusResult)) +
−
lonR <- shp2@coords[ , 1]
+
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"
 +
)
  
−
mymarkers <- data.frame(
+
########### Concentrations in time
−
lat = c(shp2@coords[,2]),
 
−
lon = c(shp2@coords[,1]),
 
−
color = colcode
 
−
)
 
  
−
#get the bounding box:
+
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"
 +
)
  
−
bb <- qbbox(lat = mymarkers[,"lat"], lon = mymarkers[,"lon"])
+
</rcode>
  
−
    if (missing(zoom))
+
==== Static GoogleMaps ====
−
        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")
 
−
    }
 
  
−
MyMap <- GetMap(
+
<rcode name='staticmapstest' graphics=1 embed=1>
−
center = c(lat.center, lon.center),
 
−
zoom = zoom,
 
−
size = size,
 
−
RETURNIMAGE = RETURNIMAGE,
 
−
GRAYSCALE = GRAYSCALE,
 
−
verbose = verbose,
 
−
...
 
−
)
 
  
−
PlotOnStaticMap(MyMap)
+
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.
  
−
PlotOnStaticMap(
+
shp <- readOGR('PG:host=localhost user=postgres dbname=spatial_db','kuopio_house')
−
MyMap,
+
proj4string(shp) <- ("+init=epsg:3067") # The map projection of NLS Finland.
−
lat = mymarkers[,"lat"],
 
−
lon = mymarkers[,"lon"],
 
−
pch = pch,
 
−
cex = cex,
 
−
col = colcode,
 
−
add = T
 
−
)
 
  
−
legend(
+
#oprint(shp@data[1:100 , ])
−
"topleft",
+
 
−
legend = names(attr(colcode, "table")),
+
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
−
title = legend_title,
+
shp2 <- spTransform(shp,epsg4326String) # Convert to longitude-latitude projection.
−
fill = attr(colcode, "palette"),
 
−
cex = 1.0,
 
−
bty = "y",
 
−
bg = "white"
 
−
)
 
−
}
 
  
−
temp <- MyRmap(shp2, maptype = "mobile", legend_title = "Ikä") # bb$lonR, bb$latR
+
MyRmap(shp2, plotvar = "ika", legend_title = "Ikä", numbins = 8, pch = 19, cex = 0.3) # Draw the map.
  
 
</rcode>
 
</rcode>
Line 343: Line 422:
  
 
* [https://gist.github.com/2292662 Plotting sp-objects to Google Maps with different CRSs]
 
* [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"