smooth_map: Create a smooth map

In mtennekes/oldtmaptools: Old Thematic Map Tools

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

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.

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

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)
    	repeats <- pmax(1, metro_eck$pop2010 %/% 1e6)
    	ids <- unlist(mapply(rep, 1:nrow(metro_eck), repeats, SIMPLIFY = FALSE))

    	metro_dots <- metro_eck[ids, ]


    	# 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 - nrow(metro_dots)*1e6)

    	# combine

    	c(st_geometry(World_pop), st_geometry(metro_dots))
    }

    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)
}

## End(Not run)

mtennekes/oldtmaptools documentation built on May 11, 2019, 8:22 p.m.