Боремся с отображением изолиний поверх карты

44 spatial/02 heatmap sample/urls.txt

@@ -0,0 +1,2 @@
+http://freakonometrics.hypotheses.org/19473
+http://blog.rtwilson.com/john-snows-cholera-data-in-more-formats/

45 GIS_tests/data/sensors_calc.csv

@@ -0,0 +1,5 @@
+"num","lon","lat"
+1,37.5774651088064,55.7621310273728
+2,37.582470748435,55.7696728117627
+3,37.5838661530659,55.763386239752
+4,37.5764183231174,55.768293454839

45 GIS_tests/data/sensors_manual.csv

@@ -0,0 +1,12 @@
+lon,lat
+37.578691607470724,55.766160765720493
+37.57990362015564,55.766504551149772
+37.580025839922193,55.765874275547077
+37.577632369493983,55.764819973580607
+37.581431367237478,55.764384492709446
+37.579109191673091,55.765158040903287
+37.58182858147876,55.766160765720493
+37.580586013852198,55.767008764295831
+37.579862880233463,55.764109449653162
+37.577856439065989,55.766567578149633
+37.576776831128157,55.76609773806296

45 GIS_tests/data/sensors_manual.gpx

@@ -0,0 +1,70 @@
+<?xml version="1.0"?>
+<gpx version="1.1" creator="GDAL 2.0.2" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ogr="http://osgeo.org/gdal" xmlns="http://www.topografix.com/GPX/1/1" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd">
+<metadata><bounds minlat="55.764109449653162" minlon="37.576776831128157" maxlat="55.767008764295831" maxlon="37.581828581478760"/></metadata>                  
+<wpt lat="55.766160765720493" lon="37.578691607470724">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.766504551149772" lon="37.57990362015564">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.765874275547077" lon="37.580025839922193">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.764819973580607" lon="37.577632369493983">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.764384492709446" lon="37.581431367237478">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.765158040903287" lon="37.579109191673091">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.766160765720493" lon="37.58182858147876">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.767008764295831" lon="37.580586013852198">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.764109449653162" lon="37.579862880233463">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.766567578149633" lon="37.577856439065989">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+<wpt lat="55.76609773806296" lon="37.576776831128157">
+  <extensions>
+    <ogr:elevation>0.0</ogr:elevation>
+    <ogr:url_name>????????? ? QGIS</ogr:url_name>
+  </extensions>
+</wpt>
+</gpx>

45 GIS_tests/data/sensors_zoo.csv

@@ -1,5 +1,6 @@
-type, lat, lon
-Center, 55.76559, 37.58022
-Sensor 1, 55.860510, 37.655286
-Sensor 2, 55.760510, 37.455286
-Sensor 3, 55.960510, 37.555286
+type, lat, lon, val
+Center, 55.76515, 37.58020, 35
+Sensor 1, 55.76364, 37.57658, 78
+Sensor 2, 55.76345, 37.58409, 46
+Sensor 3, 55.76924, 37.58356, 15
+Sensor 4, 55.76778, 37.57766, 54

45 GIS_tests/map_field.R

@@ -8,11 +8,17 @@ library(scales)
 library(readr) #Hadley Wickham, http://blog.rstudio.org/2015/04/09/readr-0-1-0/
 library(sp)       # spatial operations library
 library(leaflet)
+library(KernSmooth)
+library(akima)
+library(rgl)
 
 options(warn=1) # http://stackoverflow.com/questions/11239428/how-to-change-warn-setting-in-r
+# options(viewer = NULL) # send output directly to browser
 
 
-ifilename <- ".\\data\\sensors.csv"
+ifilename <- ".\\data\\sensors_zoo.csv"
+ifilename <- ".\\data\\sensors_calc.csv"
+ifilename <- ".\\data\\sensors_manual.csv"
 # ïîäãðóæàåì äàííûå ïî ñåíñîðàì
 mydata <- read_delim(ifilename, delim = ",", quote = "\"",
                      col_names = TRUE,
@@ -28,7 +34,8 @@ mydata <- read_delim(ifilename, delim = ",", quote = "\"",
 # for function ‘coordinates’ for signature ‘"tbl_df"’
 
 mydata <- as.data.frame(mydata)
-
+# Äîáàâèì êîëîíêó ïðèíóäèòåëüíî. Àêòóàëüíî äëÿ çàãðóçêè ðàñ÷åòíûõ ïîçèöèé ñåíñîðîâ
+mydata$val <- runif(nrow(mydata), 10, 90) # íàêèäûâàåì âëàæíîñòü â ïðîöåíòàõ
 
 # Ñíà÷àëà íàäî ïðåîáðàçîâàòü â Spatial Data
 # 1. Îïðåäåëÿåì CRS
@@ -44,12 +51,81 @@ plot(mydataSp, main = '
 # =================== leaflet ==================
 # äëÿ leaflet íåîáÿçàòåëüíà ðàáîòà ñ Spatial Objects
 #build the map without "pipes"
-m <- leaflet(mydata, base.map="water") #setup map
-m <- addTiles(m) #add open street map data
-m <- setView(m, lng = 38.6, lat = 54.9, zoom = 12)
+# m <- leaflet(mydata, base.map="osm") #setup map
+m <- leaflet(mydata) %>% #setup map
+  addTiles() %>% #add open street map data
+  setView(lng = 37.58, lat = 55.76, zoom = 15) %>%
+  addMarkers(lng = ~lon, lat = ~lat)
 # map$tileLayer(provider = 'Stamen.Watercolor') # https://github.com/ramnathv/rMaps
 #m$tileLayer(provider = 'Stamen.Watercolor')
 
 #m <- addCircleMarkers(radius = ~size, color = ~color, fill = FALSE)
+# èìåíà ïðîâàéäåðîâ ñìîòðèì â ò.÷. çäåñü: http://leaflet-extras.github.io/leaflet-providers/preview/
+# m <- addProviderTiles(m, "CartoDB.Positron")
+#m <- addProviderTiles(m, "Stamen.Watercolor")
+# m <- addProviderTiles(m, "OpenWeatherMap.RainClassic")
+
+
+m
+
+
+stop()
+# òåïåðü ïî íîâîé, áåðåì èäåè îòñþäà: http://stackoverflow.com/questions/24410292/how-to-improve-interp-with-akima
+
+fld <- with(mydata, interp(lon, lat, val,
+    xo = seq(min(lon), max(lon), length = 500),
+    yo = seq(min(lat), max(lat), length = 500)
+  )) # http://www.statmethods.net/stats/withby.html
+filled.contour(x=fld$x, y=fld$y, z=fld$z,
+               color.palette=colorRampPalette(c("white", "blue")))
+
+stop()
+
+# http://docs.ggplot2.org/0.9.3.1/stat_contour.html
+# Basic plot
+v <- ggplot(mydata, aes(x = lon, y = lat, z = val))
+# v + stat_contour(bins = 5)
+v
+
+qplot(lon, lat, z = val, data = mydata, geom = "contour")
+stop()
+
+
+s <-  interp(mydata$lon, mydata$lat, mydata$val, linear = FALSE, extrap = TRUE) #from akima
+# If linear is TRUE (default), linear interpolation is used in the triangles bounded by data points.
+# Cubic interpolation is done if linear is set to FALSE. If extrap is FALSE, z-values for points outside
+# the convex hull are returned as NA. No extrapolation can be performed for the linear case.
+mm <- s$z
+# contour(s$x, s$y, s$z)
+contour(s)
+
+CL <- contourLines(s$x, s$y, s$z, nlevels = 9)
+# nlevels óêàçûâàåò ðàçáèåíèå, íî ðåàëüíî íàðåçîê ìîæåò áûòü áîëüøå (íåñêîëüêî ãîð)
+
+m = leaflet() %>% addTiles() %>%
+  # setView(lng = 37.58, lat = 55.76, zoom = 14) %>%
+  addPolygons(CL[[1]]$x,CL[[1]]$y,fillColor = "red", stroke = FALSE) %>%
+  addPolygons(CL[[3]]$x,CL[[3]]$y,fillColor = "green", stroke = FALSE) %>%
+  addPolygons(CL[[5]]$x,CL[[5]]$y,fillColor = "blue", stroke = FALSE) %>%
+  addPolygons(CL[[7]]$x,CL[[7]]$y,fillColor = "red", stroke = FALSE) %>%
+  addPolygons(CL[[9]]$x,CL[[9]]$y,fillColor = "red", stroke = FALSE)
+
+w <- s$z
 
 m
+
+stop()
+
+# http://freakonometrics.hypotheses.org/19473
+# http://gis.stackexchange.com/questions/168886/r-how-to-build-heatmap-with-the-leaflet-package
+X <- mydata[c("lat", "lon")]
+# bw.ucv -- Bandwidth Selectors for Kernel Density Estimation
+kde2d <- bkde2D(X, bandwidth = c(bw.ucv(X[,1]), bw.ucv(X[,2])))
+# bkde2d íå î÷åíü ãîäèòñÿ, ïîñêîëüêó ýòî density estimate
+
+t <- kde2d$fhat
+
+contour(x=kde2d$x1, y=kde2d$x2, z=kde2d$fhat)
+stop()
+data(quakes)
+

45 GIS_tests/sensors.gpx

@@ -0,0 +1,47 @@
+<?xml version="1.0" encoding="utf-8"?>
+<gpx version="1.0" creator="QGIS">
+<wpt lat="55.766160765720" lon="37.578691607471">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.766504551150" lon="37.579903620156">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.765874275547" lon="37.580025839922">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.764819973581" lon="37.577632369494">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.764384492709" lon="37.581431367237">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.765158040903" lon="37.579109191673">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.766160765720" lon="37.581828581479">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.767008764296" lon="37.580586013852">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.764109449653" lon="37.579862880233">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.766567578150" lon="37.577856439066">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+<wpt lat="55.766097738063" lon="37.576776831128">
+<urlname>Оцифрован в QGIS</urlname>
+<ele>0</ele>
+</wpt>
+</gpx>