tmaptools (version 2.0-1)

smooth_map: Create a smooth map

Description

Create a smooth map from a shape object. A 2D kernel density estimator is applied to the shape, which can be a spatial points, polygons, or raster object. Various format are returned: a smooth raster, contour lines, and polygons. The covered area can be specified, i.e., the area outside of it is extracted from the output. Note that this function supports sf objects, but still uses sp-based methods (see details).

Usage

smooth_map(shp, var = NULL, nrow = NA, ncol = NA, N = 250000,
  unit = "km", unit.size = 1000, smooth.raster = TRUE, nlevels = 5,
  style = ifelse(is.null(breaks), "pretty", "fixed"), breaks = NULL,
  bandwidth = NA, threshold = 0, cover.type = NA, cover = NULL,
  cover.threshold = 0.6, weight = 1, extracting.method = "full",
  buffer.width = NA, to.Raster = NULL)

Arguments

shp

shape object of class Spatial, Raster, or sf. Spatial points, polygons, and grids are supported. Spatial lines are not.

var

variable name. Not needed for SpatialPoints. If missing, the first variable name is taken. For polygons, the variable should contain densities, not absolute numbers.

nrow

number of rows in the raster that is used to smooth the shape object. Only applicable if shp is not a SpatialGrid(DataFrame) or Raster

ncol

number of rows in the raster that is used to smooth the shape object. Only applicable if shp is not a SpatialGrid(DataFrame) or Raster

N

preferred number of points in the raster that is used to smooth the shape object. Only applicable if shp is not a SpatialGrid(DataFrame) or Raster

unit

unit specification. Needed when calculating density values. When set to NA, the densities values are based on the dimensions of the raster (defined by nrow and ncol). See also unit.size.

unit.size

size of the unit in terms of coordinate units. The coordinate system of many projections is approximately in meters while thematic maps typically range many kilometers, so by default unit="km" and unit.size=1000 (meaning 1 kilometer equals 1000 coordinate units).

smooth.raster

logical that determines whether 2D kernel density smoothing is applied to the raster shape object. Not applicable when shp is a SpatialPoints object (since it already requires a 2D kernel density estimator). Other spatial objects are converted to a raster, which is smoothed when smooth.raster=TRUE.

nlevels

preferred number of levels

style

method to cut the color scale: e.g. "fixed", "equal", "pretty", "quantile", or "kmeans". See the details in classIntervals.

breaks

in case style=="fixed", breaks should be specified

bandwidth

single numeric value or vector of two numeric values that specifiy the bandwidth of the kernal density estimator. By default, it is 1/50th of the shortest side in units (specified with unit.size).

threshold

threshold value when a 2D kernel density is applied. Density values below this threshold will be set to NA. Only applicable when shp is a SpatialPoints or smooth.raster=TRUE.

cover.type

character value that specifies the type of raster cover, in other words, how the boundaries are specified. Options: "original" uses the same boundaries as shp (default for polygons), "smooth" calculates a smooth boundary based on the 2D kernal density (determined by smooth_raster_cover), "rect" uses the bounding box of shp as boundaries (default for spatial points and grids).

cover

SpatialPolygons shape that determines the covered area in which the contour lines are placed. If specified, cover.type is ignored.

cover.threshold

numeric value between 0 and 1 that determines which part of the estimated 2D kernal density is returned as cover. Only applicable when cover.type="smooth".

weight

single number that specifies the weight of a single point. Only applicable if shp is a SpatialPoints object.

extracting.method

Method of how coordinates are extracted from the kernel density polygons. Options are: "full" (default), "grid", and "single". See details. For the slowest method "full", extract is used. For "grid", points on a grid layout are selected that intersect with the polygon. For "simple", a simple point is generated with gPointOnSurface.

buffer.width

Buffer width of the iso lines to cut kernel density polygons. Should be small enough to let the polygons touch each other without space in between. However, too low values may cause geometric errors.

to.Raster

not used anymore, since the "raster" output is always a RasterLayer as of version 2.0

Value

List with the following items:

"raster"

A smooth raster, which is either a SpatialGridDataFrame or a RasterLayer (see to.Raster)

"iso"

Contour lines, which is an sf object of spatial lines.

"polygons"

Kernel density polygons, which is an sf object of spatial polygons

"bbox"

Bounding box of the used raster

"ncol"

Number of rows in the raster

"nrow"

Number of columns in the raster

Details

For the estimation of the 2D kernal density, code is borrowed from bkde2D. This implemention is slightly different: bkde2D takes point coordinates and applies linear binning, whereas in this function, the data is already binned, with values 1 if the values of var are not missing and 0 if values of var are missing.

This function supports sf objects, but still uses sp-based methods, from the packages sp, rgeos, and/or rgdal.

Examples

# NOT RUN {
if (require(tmap)) {
    # set mode to plotting mode
    current.mode <- tmap_mode("plot")

    ####################################
    ## Already smoothed raster
    ####################################
    vol <- raster::raster(t(volcano[, ncol(volcano):1]), xmn=0, xmx=870, ymn=0, ymx=610)
    vol_smooth <- smooth_map(vol, smooth.raster = FALSE, nlevels = 10)

    tm_shape(vol_smooth$polygons) +
    	tm_fill("level", palette=terrain.colors(11), title="Elevation") +
    	tm_shape(vol_smooth$iso) +
    	tm_iso(col = "black", size = .7, fontcolor="black") +
    	tm_layout("Maunga Whau volcano (Auckland)", title.position=c("left", "bottom"),
    	    inner.margins=0) +
    	tm_legend(width=.13, position=c("right", "top"), bg.color="gray80", frame = TRUE)


    ####################################
    ## Smooth polygons
    ####################################
    data(NLD_muni)

    NLD_muni$population_dens <- as.vector(calc_densities(NLD_muni, "population"))

    qtm(NLD_muni, fill="population_dens")

    NLD_smooth <- smooth_map(NLD_muni, var = "population_dens")

    qtm(NLD_smooth$raster, style="grey")
    qtm(NLD_smooth$polygons, format="NLD")

    ####################################
    ## Smooth points
    ####################################

    # Approximate world population density as spatial points, one for each 1 million people,
    # in the following way. Each metropolitan area of x million people will be represented
    # by x dots. The remaining population per country will be represented by dots that are
    # sampled across the country.
    create_dot_per_1mln_people <- function() {
    	data(World, metro)
    	metro_eck <- set_projection(metro, projection = "eck4")

    	# aggregate metropolitan population per country
    	metro_per_country <- tapply(metro_eck$pop2010, INDEX = list(metro_eck$iso_a3), FUN=sum)
    	metro_per_country_in_World <- metro_per_country[names(metro_per_country) %in% World$iso_a3]

    	# assign to World shape
    	World$pop_metro <- 0
    	World$pop_metro[match(names(metro_per_country_in_World), World$iso_a3)] <-
    		metro_per_country_in_World

    	# define population density other than metropolitan areas
    	World$pop_est_dens_non_metro <- (World$pop_est - World$pop_metro) / World$area

    	# generate dots for metropolitan areas (1 dot = 1mln people)
    	metro_dots <- do.call("sbind", lapply(1:length(metro_eck), function(i) {
    		m <- metro_eck[i,]
    		m[rep(1, max(1, m$pop2010 %/% 1e6)),]
    	}))

    	# sample population dots from non-metropolitan areas (1 dot = 1mln people)
    	World_pop <- sample_dots(World, vars="pop_est_dens_non_metro", w = 1e6,
    							 npop = 7.3e9 - length(metro_dots)*1e6)

    	# combine
    	sbind(as(World_pop, "SpatialPoints"), as(metro_dots, "SpatialPoints"))
    }

    World_1mln_dots <- create_dot_per_1mln_people()


    # dot map
    tm_shape(World_1mln_dots) + tm_dots()

    # create smooth map
    World_list <- smooth_map(World_1mln_dots, cover = World, weight=1e6)

    # plot smooth raster map
    qtm(World_list$raster, style="grey")

    # plot smooth raster map
    qtm(World, bbox="India") + qtm(World_list$iso)

    # plot kernel density map
    qtm(World_list$polygons, style="grey", format="World")

    ####################################
    ## Smooth raster
    ####################################
    data(land)

    land_smooth <- smooth_map(land, var="trees", cover.type = "smooth")

    qtm(land, raster="trees")
    qtm(land_smooth$raster)
    qtm(land_smooth$polygons, format="World", style="grey")

    # reset current mode
    tmap_mode(current.mode)
}
# }
# NOT RUN {
# }