Opasnet map

[show] 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 }}

{{#display_map: 62.928986, 28.706245; 62.937734, 28.712769 }}

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.

See also

• Plotting sp-objects to Google Maps with different CRSs
• Map editor
• Batch geocoding: give addresses, the website will provide geolocations. Works for a few hundred addresses.

Related files

<mfanonymousfilelist></mfanonymousfilelist>