If you are not redirected automatically, follow the link to new page mapview - interactive viewing of spatial data in R

(CC-BY-NC-SA)

This document was produced on Di Dez 15 2015 using mapview version 1.0.0

Installation

The easiest way to install mapview is to use library("devtools"). With devtools installed and loaded type the following:

For the stable release of mapview:

install_github("environmentalinformatics-marburg/mapview")

For the devolpment version:

install_github("environmentalinformatics-marburg/mapview", ref = "develop")

Introduction

Working with spatial data in R I find myself quite often in the need to quickly visually check whether a certain analysis has produced reasonable results. There are two ways I usually do this. Either I:

  1. (sp)plot the data in R and then toggle back and forth between the static plots (I use RStudio) or
  2. save the data to the disk and then open in QGIS or similar to interactively examine the results.

Both these approaches are semi-optimal. Where option 1. is fine for a quick glance at a coarse patterns, it lacks the possibility to have a closer look into the results via zooming and panning. While option 2. provides the interactivity, the detour via the hard disk is annoying (at best), especially when fine-tuning and checking regularly.

I attended this years useR2015! conference in Aalborg (which was marvelous!) and attended the session on interactive graphics in R where Joe Cheng from RStudio presented the leaflet package. leaflet is great as it gives you a great deal of control over the individual map components. This, however, also means that in order to get some useful visualization of spatial objects, we need to do a fair bit of coding for the map components to show what we want (e.g. we would need to provide popup-text manually for each object that we want to map). What an interactive functionality would need is some default behaviour for different objects from the spatial universe.

This got me thinking and sparked my enthusiasm to write some wrapper functions for leaflet to provide at least very basic interactive graphing capabilities that are directly accessible within RStudio (or the web browser, if you’re not using RStudio).

The result is a package called mapview. The main workhorse function is mapView() (mind the capital ‘V’) and is currently defined for:

General design

A call to mapView() will return an object of class mapview. This class has 2 slots

  • @object - a list of the objects that are displayed on the map. This means that this slot will contain the re-projected (and in the case of Raster objects possibly re-sampled) objects which enables tracing of the modifications that took place.
  • @map - the leaflet map. This is an S3 class object (see leaflet documentation for details on the specifics).

By default mapView() provides two base layers between which one can toggle (a standard OpenStreeMap layer and ESRI’s WorldImagery layer). Depending on the object and/or argument settings one or several layers are created (each with or without it’s own legend) that can also be toggled. Note that in order to render properly, all layers need to be re-projected to leaflet’s underlying web mercator projection

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs

which in the case of large Raster* objects or Spatial objects with lots of features can be time consuming.

In the following I would like to present a few use case scenarios that highlight the current capabilities of mapview (largely taken from the help files):

Raster data objects

NOTE: similar to spplot() for Raster* objects mapView() has a maxpixels = argument to avoid long rendering times for large rasters. This is currently set to a default value of 500000 (which produces acceptable rendering times on my machine) and a suitable warning is printed for larger Raster objects.

RasterLayer

library(mapview)
library(raster)

mapView(poppendorf[[10]])

RasterStack/Brick

If you pass a RasterStack/Brick to mapView() this will create one map layer + legend for each RasterLayer. This is likely to cause problems for larger Stacks as there is currently no way to tell mapView() to only show the legend of the highlighted layer. This means that some of the legends won’t show in the viewer.

mapView(poppendorf)

There is (to my knowledge) currently no way to provide this functionality in the leaflet package, though I have seen solutions for this in JavaScript which means that this will likely be available in the future.

SpatialPixelsDataFrame

Given that a SpatialPixelsDataFrame has an attribute table we can use argument zcol = "column_name" to plot specific columns of the table. The zcol argument can be used with the …DataFrame versions of all Spatial* objects that are currently supported (see below).

library(sp)

data(meuse.grid)
coordinates(meuse.grid) <- ~x+y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
gridded(meuse.grid) <- TRUE

mapView(meuse.grid, zcol = "soil")

Vector data objects

For spatial vector objects from the sp-package I would like to highlight the following arguments:

SpatialPoints(DataFrame)

data(meuse)
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")

# only one layer, all info in popups
mapView(meuse)

Another argument you can use with SpatialPointsDataFrame is cex which takes either an integer value (default is 10) or the character name of one of the columns in the attribute table. If the latter is provided, circles are scaled relative to the respective values in the column provided.

mapView(meuse, cex = "lead")

Apart from the cex argument, the same functionality is implemented for

SpatialPolygons(DataFrame)

data("gadmCHE")
mapView(gadmCHE)

and

SpatialLines(DataFrame)

data("atlStorms2005")
mapView(atlStorms2005, burst = TRUE)

Additional functionality

For the remote sensing inclined, mapview-package offers some additional functionality:

Red-Green-Blue images

viewRGB() provides the possibility to view Red-Green-Blue true-/false-color images of RasterStacks/Bricks. Any band combination is possible.

True color example

data("poppendorf")
viewRGB(poppendorf, 4, 3, 2) #Bands 4,3,2 correspond to Red Green Blue

Note, by default the color to render NA values is set to "transparent", so let’s change that to "black"

False color example

viewRGB(poppendorf, 5, 4, 3, na.color = "black") #Bands 5,4,3 correspond to NIR Red Green

Fun fact: within this tiny area one can find at least five breweries (see below)… Welcome to the land of Beer!

View extent only

As mentioned above, it is possible to view only the extent of Raster and bbox of Spatial objects. Clicking on the rectangle triggers a pop-up window showing the coordinates of the corners of the extent in epsg:4326.

Raster* object extent

viewExtent(poppendorf)

Spatial* object bbox

viewExtent(meuse, color = "black")

Overlays

To make mapview even more user-friendly, I have defined a +-method for combinations of

  • mapview + mapview
  • mapview + ANY (ANY here refers to all classes supported by mapView())
  • leaflet + ANY (so that one can easily pass spatial objects to existing leaflet maps)

This means you can easily combine objects by enabling things like

mapView(meuse.grid, zcol = "soil") + meuse

Basically, any combination is possible.

The zoom will be set to the extent of the objected that was added last. Also, layers are overlaid according to the order of calls. Note, however, that leaflet will always render vector data on top of raster data,

data("breweries91")
mapView(breweries91, color = "red") + viewRGB(poppendorf, 4, 3, 2) + viewExtent(poppendorf)

So here you have it, 5 breweries within the extent of the RGB image

Non-projected data

It is possible to view objects that are not projected (i.e. have their coordinate reference system - CRS - set to NA). If this is the case, coordinates are rescaled to lie between 0 and 1 and then the y-coordinates are multiplied by the ratio of y-extent / x-extent. This means that x-coordinates will always have min == 0 and max == 1 while y-coordinates will have min == 0 but their max will depend on the ratio, so that this can be slightly higher or lower than 1.

When we supply non-projected data to mapView() a suitable warning will be printed to the console (see below).

Here’s an example for point data. For this we load the meuse data set once again, but this time we won’t set the CRS.

data("meuse")
coordinates(meuse) <- ~x+y

## check coordinate reference system
proj4string(meuse)
## [1] NA

Yet, even though it has no CRS, we can still visulize it. Though we obviously cannot provide any background maps for the visualization.

mapView(meuse)
## Warning in spCheckAdjustProjection(x): supplied SpatialPointsDataFrame has no projection information!
##  scaling coordinates and showing layer without background map

Note, how the coordinates have been rescaled.

This is also possible for Raster* objects, which allows us to view pretty much any type of raster data, including those that cannot be projected, such as photos. This is quite handy as it enables closer investigation of non-standard spatial data, for example gathered through imaging ground-based remote sensing techniques (e.g. thermal imaging).

Here’s a toy example using one of the header images from the Environmental Informtics Marburg web page that shows the author of this package during an R teaching session at our research station at Mt. Kilimanjaro in Tanzania.

## get image from the web page
library(jpeg)

web_image <- "http://umweltinformatik-marburg.de/uploads/tx_rzslider/teaching_header_kili_resized.jpg"

jpg <- readJPEG(readBin(web_image, "raw", 1e6))

## Convert image data to raster
rst_blue <- raster(jpg[, , 1])
rst_green <- raster(jpg[, , 2])
rst_red <- raster(jpg[, , 3])

img <- brick(rst_red, rst_green, rst_blue)

## visualize image in viewer
viewRGB(img)
## Warning in rasterCheckAdjustProjection(x): supplied RasterBrick has no projection information!
##  scaling coordinates and showing layer without background map

In case you are interested, under the hood, the rendering of unprojected Raster* data works similar to that of Spatial* objects. Corner coordinates of the extent are rescaled (y-coordinates are multiplied with the ratio of nrow / ncol to preserve dimensions). In this case, however, we need to provide a CRS as leaflet expects one. Therefore, we set this to standard latitude/longitude EPSG:4236. This means, that the rendered image actually does have a spatial reference but mapview simply does not provide basemaps. If we were to provide one, we would see that our data lies somewhere in the Gulf of Guinea.

Applications for non-projected data are obviously limited, yet, especially the possibility to view standard images interactively provides a nice tool for checking parts of an image(stack) in more detail.

Keep in mind that speed and memory limitations for Raster* data still applies, though speed will be much less of an issue, as we do not need to reporject the data.

leaflet integration

Integrating mapview and leaflet works both ways…

leaflet + mapview

You can first create a map using leaflet and then add spatial objects as you like with mapview allowing you to flexibly create maps in whatever way you like.

## first set the CRS again
proj4string(meuse) <- CRS("+init=epsg:28992")

m <- leaflet() %>%
  addProviderTiles("Stamen.Toner")
m + meuse

mapview + leaflet

Given that objects created with mapView() contain a slot called map, you can first create a map using mapView() and then add components provided by leaflet. To reproduce the following example you will (at the time of this writing) need to install leaflet from the GitHub leaflet repository which provides new functionality to measure distances and areas. The following code will then add a control element in the top right corner of the map which lets you interactively create a polygon and will show its perimeter and area. Double-click to finish a measurement.

m <- mapView(breweries91) + viewExtent(poppendorf)
m@map %>% addMeasure(primaryLengthUnit = "meters", primaryAreaUnit = "sqmeters",
                     activeColor = "#3D535D", completedColor = "#7D4479")

So now you can measure the distances between the breweries or figure out how big the area is that is covered by the extent…

Advanced example

As a final, more advanced, example of what can be done with mapview we will mimic a common workflow when working with spatial data. We will calculate the mean January temperature for each district of Switzerland using raster::extract. Then we will render this on top of the mean January temperatures and also add a hill-shade raster of Switzerland. All data sets are included in mapview. They were created following the examples shown in the help pages of raster::getData and raster::hillShade. This example also highlights a few arguments that can be adjusted which were not mentioned so far.

data("tmeanCHE")
data("hillshadeCHE")
data("gadmCHE")

t <- extract(tmeanCHE, gadmCHE, fun = mean, na.rm = TRUE, sp = TRUE)
t$Jan.mean.T <- round(t$Jan.mean.T, 2)

m1 <- mapView(hillshadeCHE, col.regions = grey.colors(10000),
              alpha.regions = 1, layer.name = "Hillshade",
              map.types = "Acetate.hillshading", legend = FALSE)
m2 <- mapView(tmeanCHE, alpha.regions = 0.6,
              map.types = c("OpenStreetMap.BlackAndWhite",
                            "OpenStreetMap.HOT"),
              layer.name = "JanTemp")
m3 <- mapView(t, zcol = "Jan.mean.T", alpha.regions = 0.9)
m1 + m2 + m3

Note how we

I hope this highlights that despite providing sensible default behaviour for viewing spatial data mapview also provides a fair deal of flexibility.

And for those who…

… still wonder where exactly the difference between mapview and leaflet lies, consider the following:

Let’s reproduce a simple mapView() call using leaflet()

Here’s the mapview version:

mapView(meuse)

And now let’s get the identical result using leaflet:

## first we need to reproject meuse to geographic coordinates
meuse <- spTransform(meuse,
                     CRSobj = CRS("+proj=longlat +datum=WGS84 +no_defs"))

## then we need to generate the text for the pop-ups
df <- as.data.frame(sapply(meuse@data, as.character),
                    stringsAsFactors = FALSE)

nms <- names(df)

txt_x <- paste0("x: ", round(coordinates(meuse)[, 1], 2))
txt_y <- paste0("y: ", round(coordinates(meuse)[, 2], 2))

txt <- rbind(sapply(seq(nrow(meuse@data)), function(i) {
  paste(nms, df[i, ], sep = ": ")
}), txt_x, txt_y)

txt <- sapply(seq(ncol(txt)), function(j) {
  paste(txt[, j], collapse = " <br/> ")
})

## finally we can create our map
leaflet(meuse) %>%
  addProviderTiles("OpenStreetMap", group = "OpenStreetMap") %>%
  addProviderTiles("Esri.WorldImagery", group = "Esri.WorldImagery") %>%
  addCircleMarkers(lng = coordinates(meuse)[, 1],
                   lat = coordinates(meuse)[, 2],
                   group = "meuse",
                   color = mapViewPalette(3)[length(mapViewPalette(3))],
                   popup = txt) %>%
  addLayersControl(position = "bottomleft",
                   baseGroups = c("OpenStreetMap",
                                 "Esri.WorldImagery"),
                   overlayGroups = "meuse")

I guess this highlights quite well what mapview is intended for. It serves to save you some effort to interactively examine your spatial data by giving you the possibility to render it using a one-liner of code.

Final thoughts

Where to go from here?

Well, there are still a lot of limitations that need addressing:

In future releases I would like to

I hope that mapview will prove useful for some people out there.

If you have any feedback, please don’t hesitate to contact me.

Bug reports should be filed at https://github.com/environmentalinformatics-marburg/mapview/issues

Best,

Tim

?>