A step beyond creating kml files of your digital soil information is the creation of customized interactive mapping products that can be visualized within your web browser. Some notes about how to do this using R are given in our book and with example here. I wanted to teach myself how i might be able to embed these interactive maps into a website, so i created an example to help me. In doing this i hope others who have not come across interactive mapping before are able to get some basic insights about it.
Interactive mapping makes sharing your data with colleagues simpler, and importantly improves the visualization experience via customization features that are difficult to achieve via the Google Earth software platform. The interactive mapping is made possible via the Leaflet R package. Leaflet is one of the most popular open-source JavaScript libraries for interactive maps. The Leaflet R package makes it easy to integrate and control Leaflet maps in R. More detailed information about Leaflet can be found here, and information specifically about the R package is here.
There is a common workflow for creating Leaflet maps in R. First is the
creation of a map widget (calling ); followed by the adding of layers or
features to the map by using layer functions (e.g. addTiles,
addMarkers, addPolygons) to modify the map widget. The map can then
be printed and visualized in the R image window or saved to HTML file
for visualization within a web browser or even embed it into a website
like the one you are reading. The following R script is a quick taste of
creating an interactive Leaflet map. It is assumed that the leaflet
and magrittr packages are installed.
First lets get our workflow initialised by loading up our R packages.
# Required R Packages
library(leaflet)
library(magrittr)
library(ithir)
library(sp)
library(raster)
library(RColorBrewer)
library(rgdal)
Point Data
We will be working with a small data set of soil information that was collected from the Hunter Valley, NSW in 2010 called . This data set is contained in the package. So first load it in:
data(HV100)
str(HV100)
## 'data.frame': 100 obs. of 6 variables:
## $ site: Factor w/ 100 levels "a1","a10","a11",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ x : num 337860 344060 347035 338235 341760 ...
## $ y : num 6372416 6376716 6372741 6368766 6366016 ...
## $ OC : num 2.03 2.6 3.42 4.1 3.04 4.07 2.95 3.1 4.59 1.77 ...
## $ EC : num 0.129 0.085 0.036 0.081 0.104 0.138 0.07 0.097 0.114 0.031 ...
## $ pH : num 6.9 5.1 5.9 6.3 6.1 6.4 5.9 5.5 5.7 6 ...
Using the coordinates function from the sp package we can define
which columns in the data frame refer to actual spatial coordinates—
here the coordinates are listed in columns x and y.
coordinates(HV100) <- ~x + y
str(HV100)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 100 obs. of 4 variables:
## .. ..$ site: Factor w/ 100 levels "a1","a10","a11",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ OC : num [1:100] 2.03 2.6 3.42 4.1 3.04 4.07 2.95 3.1 4.59 1.77 ...
## .. ..$ EC : num [1:100] 0.129 0.085 0.036 0.081 0.104 0.138 0.07 0.097 0.114 0.031 ...
## .. ..$ pH : num [1:100] 6.9 5.1 5.9 6.3 6.1 6.4 5.9 5.5 5.7 6 ...
## ..@ coords.nrs : int [1:2] 2 3
## ..@ coords : num [1:100, 1:2] 337860 344060 347035 338235 341760 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:100] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "x" "y"
## ..@ bbox : num [1:2, 1:2] 335160 6365091 350960 6382816
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NA
Note now that by using the str function, the class of HV100 has now
changed from a dataframe to a SpatialPointsDataFrame. The
SpatialPointsDataFrame structure is essentially the same data frame,
except that additional spatial elements have been added or partitioned into slots. Some important ones being the bounding box (sort of like the spatial extent of the data), and the coordinate reference system (proj4string), which we need to define for our data set. To define the CRS, we have to know some things about where our data are from, and what was the corresponding CRS used when recording the spatial information in the field. For this data set the CRS used was WGS1984 UTM Zone 56. To explicitly tell R this information we define the CRS as a character string which describes a reference system in a way understood by the PROJ.4 projection library. An interface to the PROJ.4 library is available in thergdalpackage. Alternative to using Proj4 character strings, we can use the corresponding yet simpler EPSG code (European Petroleum Survey Group).rgdal
also recognizes these codes. If you are unsure of the Proj4 or EPSG code
for the spatial data that you have, but know the CRS, you should consult
Spatial Reference for assistance. The
EPSG code for WGS1984 UTM Zone 56 is: 32556. So lets define to CRS for
this data.
proj4string(HV100) <- CRS("+init=epsg:32756")
HV100@proj4string
## CRS arguments:
## +init=epsg:32756 +proj=utm +zone=56 +south +datum=WGS84 +units=m
## +no_defs +ellps=WGS84 +towgs84=0,0,0
To look at the locations of the data in Google Earth, we first need to
make sure the data is in the WGS84 geographic CRS. If the data is not in
this CRS (which is not the case for this data), then we need to perform
a coordinate transformation. This is facilitated by using the
spTransform function in sp. The EPSG code for WGS84 geographic is:
- We can then export out our transformed HV100 data set to a KML file and visualize it in Google Earth. We will use it later for our interactive web mapping.
HV100.ll <- spTransform(HV100, CRS("+init=epsg:4326"))
#Export KML
#writeOGR(HV100.ll, "HV100.kml", "ID", "KML")
Rasters
Most of the functions needed for handling raster data are contained in
the raster package. There are functions for reading and writing raster
files from and to different raster formats. In DSM we work quite a deal
with data in table format and then rasterise this data so that we can
make a map. To do this in R, lets bring in a data frame. This could be
either from a text-file, but as for the previous occasions the data is
imported from the ithir package. This data is a digital elevation
model with 100m grid resolution, from the Hunter Valley, NSW, Australia.
data(HV_dem)
str(HV_dem)
## 'data.frame': 21518 obs. of 3 variables:
## $ X : num 340210 340310 340110 340210 340310 ...
## $ Y : num 6362641 6362641 6362741 6362741 6362741 ...
## $ elevation: num 177 175 178 172 173 ...
As the data is already a raster (such that the row observation indicate
locations on a regular spaced grid), but in a table format, we can just
use the rasterFromXYZ function from raster. Also we can define the
CRS just like we did with the HV100 point data we worked with before.
r.DEM <- rasterFromXYZ(HV_dem)
proj4string(r.DEM) <- CRS("+init=epsg:32756")
For web mapping and associated Google Earth mapping, remember that we
need to reproject our data because it is in the UTM system, and need to
get it to WGS84 geographic. The raster re-projection is performed using
the projectRaster function. Look at the help file for this function.
Probably the most important parameters are crs, which takes the CRS
string of the projection you want to convert the existing raster to,
assuming it already has a defined CRS. The other is method which
controls the interpolation method. For continuous data,
bilinear would be suitable, but for categorical,ngb, (which is
nearest neighbor interpolation) is probably better suited.
p.r.DEM <- projectRaster(r.DEM, crs = "+init=epsg:4326", method = "bilinear")
Interactive Mapping
Now we have our data. Lets first provide an example of what is meant by an interactive map. Explaination is to follow.
leaflet() %>%
addMarkers(lng = 151.210558, lat = -33.852543,
popup = "The view from here is amazing!") %>%
addProviderTiles("Esri.WorldImagery")
Interactive features of this map include markers with text, plus ability
to zoom and map panning. More will be discussed about the layer
functions of the leaflet map further on. What has not been encountered
yet is the forward pipe operator %>%. This operator will forward a
value, or the result of an expression, into the next function call or
expression. To use this operator the magrittr package is required. The
example script below shows the same example using and not using the
forward pipe operator.
#Draw 100 random uniformly distributed numbers between 0 and 1
x <- runif(100)
sqrt(sum(range(x)))
## [1] 1.005229
##..is equivalent to (using forward pipe operator)
x %>% range %>% sum %>% sqrt
## [1] 1.005229
Sometimes what we want to do in R can get lost within a jumble of
brackets, whereas using the forward pipe operator the process of
operations is a lot clearer. So lets begin to construct some Leaflet
mapping using our prepared data from a little earlier regarding the
point (HV100.ll) and raster data (p.r.DEM). Firstly, lets create a
basic map — example of not using and then using the forward pipe
operator.
# Basic map
#without piping operator
addMarkers(addTiles(leaflet()), data = HV100.ll)
#with forward pipe operator
leaflet() %>%
addTiles() %>%
addMarkers(data = HV100.ll)
With the above, we are calling upon a pre-existing base map via the
addTiles() function. Leaflet supports base maps using map tiles,
popularized by Google Maps and now used by nearly all interactive web
maps. By default,
OpenStreetMap
tiles are used. Alternatively, many popular free third-party base maps
can be added using the addProviderTiles() function, which is
implemented using the leaflet-providers plugin. For example, previously
we used the Esri.WorldImagery base mapping. The full set of possible
base maps can be found
here.
Note that an internet connection is required for access to the base maps
and map tiling. The last function used above the the addMarkers
function, we we simply call up the point data we used previously, which
are those soil point observations and measurements from the Hunter
Valley, NSW. A basic map will have been created with your plot window.
For the next step, lets populate the markers we have created with some
of the data that was measured, then add the Esri.WorldImagery base
mapping.
# Populate pop-ups
my_pops <- paste0(
"<strong>Site: </strong>",
HV100.ll$site,
'<br>
<strong> Organic Carbon (%): </strong>',
HV100.ll$OC,
'<br>
<strong> soil pH: </strong>',
HV100.ll$pH
)
# Create interactive map
leaflet() %>%
addProviderTiles("Esri.WorldImagery") %>%
addMarkers(data = HV100.ll, popup = my_pops)
Further, we can colour the markers and add a map legend. Here we will
get the quantiles of the measured SOC percentage and color the markers
accordingly. Note that you will need the colour ramp package
RColorBrewer installed.
# Colour ramp
pal1 <- colorQuantile("YlOrBr", domain = HV100.ll$OC)
# Create interactive map
leaflet() %>%
addProviderTiles("Esri.WorldImagery") %>%
addCircleMarkers(data = HV100.ll, color = ~pal1(OC), popup = my_pops) %>%
addLegend("bottomright", pal = pal1, values = HV100.ll$OC,
title = "Soil Organic Carbon (%) quantiles",
opacity = 0.8
)
It is very worth consulting the help files associated with the leaflet
R package for further tips on creating further customized maps. The
website dedicated to that package, which was mentioned above is also a
very helpful resource too.
Raster maps can also be featured in our interactive mapping too, as illustrated in the following script.
#Colour ramp
pal2 <- colorNumeric(
brewer.pal(n = 9, name = "YlOrBr"),
domain = values(p.r.DEM),
na.color = "transparent"
)
#interactive map
leaflet() %>%
addProviderTiles("Esri.WorldImagery") %>%
addRasterImage(p.r.DEM, colors = pal2, opacity = 0.7) %>%
addLegend("topright", opacity=0.8, pal = pal2, values = values(p.r.DEM),
title = "Elevation")
)