Thematic Map Tools

Description

This package offers a set of handy tool functions for reading and processing spatial data. The aim of these functions is to supply the workflow to create thematic maps, e.g. read shape files, set map projections, append data, calculate areas and distances, and query OpenStreetMap. The visualization of thematic maps can be done with the tmap package.

Details

This page provides a brief overview of all package functions.

Tool functions (shape)

Tool functions (colors)

Spatial transformation functions

Input and output functions

Author(s)

Maintainer: Martijn Tennekes mtennekes@gmail.com

See Also

Useful links:

• https://github.com/r-tmap/tmaptools

• https://r-tmap.github.io/tmaptools/

• Report bugs at https://github.com/r-tmap/tmaptools/issues

Approximate area sizes of the shapes

Description

Approximate the area sizes of the polygons in real-world area units (such as sq km or sq mi), proportional numbers, or normalized numbers. Also, the areas can be calibrated to a prespecified area total. This function is a convenient wrapper around st_area.

Usage

approx_areas(shp, target = "metric", total.area = NULL)

Arguments

Details

Note that the method of determining areas is an approximation, since it depends on the used projection and the level of detail of the shape object. Projections with equal-area property are highly recommended. See https://en.wikipedia.org/wiki/List_of_map_projections for equal area world map projections.

Value

Numeric vector of area sizes (class units).

See Also

approx_distances

Examples

if (require(tmap)) {
    data(NLD_muni)

    NLD_muni$area <- approx_areas(NLD_muni, total.area = 33893)

    tm_shape(NLD_muni) +
        tm_bubbles(size="area",
                   size.legend = tm_legend(title = expression("Area in " * km^2)))


    # function that returns min, max, mean and sum of area values
    summary_areas <- function(x) {
        list(min_area=min(x),
             max_area=max(x),
             mean_area=mean(x),
             sum_area=sum(x))
    }

    # area of the polygons
    summary_areas(approx_areas(NLD_muni))

    # area of the polygons, adjusted corrected for a specified total area size
    summary_areas(approx_areas(NLD_muni, total.area=33893))

    # proportional area of the polygons
    summary_areas(approx_areas(NLD_muni, target = "prop"))

    # area in squared miles
    summary_areas(approx_areas(NLD_muni, target = "mi mi"))

    # area of the polygons when unprojected
    summary_areas(approx_areas(sf::st_transform(NLD_muni, crs = 4326)))
}

Approximate distances

Description

Approximate distances between two points or across the horizontal and vertical centerlines of a bounding box.

Usage

approx_distances(x, y = NULL, projection = NULL, target = NULL)

Arguments

Value

If y is specifyed, a list of two: unit and dist. Else, a list of three: unit, hdist (horizontal distance) and vdist (vertical distance).

See Also

approx_areas

Examples

## Not run: 
if (require(tmap)) {
    data(NLD_prov)

    # North-South and East-West distances of the Netherlands
    approx_distances(NLD_prov)

    # Distance between Maastricht and Groningen
    p_maastricht <- geocode_OSM("Maastricht")$coords
    p_groningen <- geocode_OSM("Groningen")$coords
    approx_distances(p_maastricht, p_groningen, projection = 4326, target = "km")

    # Check distances in several projections
    sapply(c(3035, 28992, 4326), function(projection) {
        p_maastricht <- geocode_OSM("Maastricht", projection = projection)$coords
        p_groningen <- geocode_OSM("Groningen", projection = projection)$coords
        approx_distances(p_maastricht, p_groningen, projection = projection)
    })
}

## End(Not run)

Bounding box generator

Description

Swiss army knife for bounding boxes. Modify an existing bounding box or create a new bounding box from scratch. See details.

Usage

bb(
  x = NA,
  ext = NULL,
  cx = NULL,
  cy = NULL,
  width = NULL,
  height = NULL,
  xlim = NULL,
  ylim = NULL,
  relative = FALSE,
  asp.target = NULL,
  asp.limit = NULL,
  current.projection = NULL,
  projection = NULL,
  output = c("bbox", "matrix", "extent")
)

Arguments

Details

An existing bounding box (defined by x) can be modified as follows:

• Using the extension factor ext.

• Changing the width and height with width and height. The argument relavitve determines whether relative or absolute values are used.

• Setting the x and y limits. The argument relavitve determines whether relative or absolute values are used.

A new bounding box can be created from scratch as follows:

• Using the extension factor ext.

• Setting the center coorinates cx and cy, together with the width and height.

• Setting the x and y limits xlim and ylim

Value

bounding box (see argument output)

See Also

geocode_OSM

Examples

if (require(tmap)) {

    ## load shapes
    data(NLD_muni)
    data(World)

    ## get bounding box (similar to sp's function bbox)
    bb(NLD_muni)

    ## extent it by factor 1.10
    bb(NLD_muni, ext=1.10)

    ## double the width
    bb(NLD_muni, width=2, relative = TRUE)

    ## crop both dimensions from 0.25 to 0.75
    bb(NLD_muni, xlim=c(.25, .75), ylim=c(.25, .75), relative = TRUE)

    ## extent it such that aspect ratio is 1
    bb(NLD_muni, asp.target = 1)

    ## convert to longlat (EPSG 4326)
    bb(NLD_muni, projection=4326)
}

## Not run: 
if (require(tmap)) {
    bb("Limburg", projection = 28992)
    bb_italy <- bb("Italy", projection = "+proj=eck4")

    tm_shape(World, bbox=bb_italy) + tm_polygons()
    # shorter alternative: tm_shape(World, bbox="Italy") + tm_polygons()
}
## End(Not run)

Convert bounding box to a spatial polygon

Description

Convert bounding box to a spatial (sfc) object . Useful for plotting (see example). The function bb_earth returns a spatial polygon of the 'boundaries' of the earth, which can also be done in other projections (if a feasible solution exists).

Usage

bb_poly(x, steps = 100, stepsize = NA, projection = NULL)

bb_earth(
  projection = NULL,
  stepsize = 1,
  earth.datum = 4326,
  bbx = c(-180, -90, 180, 90),
  buffer = 1e-06
)

Arguments

Value

sfc object

Examples

if (require(tmap) && packageVersion("tmap") >= "2.0") {
    data(NLD_muni)

    current.mode <- tmap_mode("view")
    qtm(bb_poly(NLD_muni))

    # restore mode
    tmap_mode(current.mode)
}

Calculate densities

Description

Transpose quantitative variables to densitiy variables, which are often needed for choroplets. For example, the colors of a population density map should correspond population density counts rather than absolute population numbers.

Usage

calc_densities(
  shp,
  var,
  target = "metric",
  total.area = NULL,
  suffix = NA,
  drop = TRUE
)

Arguments

Value

Vector or data.frame (depending on whether length(var)==1 with density values.

Examples

if (require(tmap) && packageVersion("tmap") >= "3.99") {
data(NLD_muni)

NLD_muni_pop_per_km2 <- calc_densities(NLD_muni,
  target = "km km", var = c("population", "dwelling_total"))
NLD_muni <- sf::st_sf(data.frame(NLD_muni, NLD_muni_pop_per_km2))

tm_shape(NLD_muni) +
  tm_polygons(
    fill = c("population_km.2", "dwelling_total_km.2"),
    fill.legend =
      list(
        tm_legend(expression("Population per " * km^2)),
        tm_legend(expression("Dwellings per " * km^2)))) +
tm_facets(free.scales = TRUE) +
  tm_layout(panel.show = FALSE)
}

Crop shape object

Description

Crop a shape object (from class sf, stars, sp, or raster). A shape file x is cropped, either by the bounding box of another shape y, or by y itself if it is a SpatialPolygons object and polygon = TRUE.

Usage

crop_shape(x, y, polygon = FALSE, ...)

Arguments

Details

This function is similar to crop from the raster package. The main difference is that crop_shape also allows to crop using a polygon instead of a rectangle.

Value

cropped shape, in the same class as x

See Also

bb

Examples

if (require(tmap) && packageVersion("tmap") >= "3.99") {
    data(World, NLD_muni, land, metro)

    #land_NLD <- crop_shape(land, NLD_muni)

    #qtm(land_NLD, raster="trees", style="natural")

    metro_Europe <- crop_shape(metro, World[World$continent == "Europe", ], polygon = TRUE)

    qtm(World) +
    tm_shape(metro_Europe) +
    	tm_bubbles("pop2010",
    	           col="red",
    	           size.legend = tm_legend("European cities")) +
    	tm_legend(frame=TRUE)
}

Geocodes a location using OpenStreetMap Nominatim

Description

Geocodes a location (based on a search query) to coordinates and a bounding box. Similar to geocode from the ggmap package. It uses OpenStreetMap Nominatim. For processing large amount of queries, please read the usage policy (https://operations.osmfoundation.org/policies/nominatim/).

Usage

geocode_OSM(
  q,
  projection = NULL,
  return.first.only = TRUE,
  keep.unfound = FALSE,
  details = FALSE,
  as.data.frame = NA,
  as.sf = FALSE,
  geometry = c("point", "bbox"),
  server = "https://nominatim.openstreetmap.org"
)

Arguments

Value

If as.sf then a sf object is returned. Else, if as.data.frame, then a data.frame is returned, else a list.

See Also

rev_geocode_OSM, bb

Examples

## Not run: 
if (require(tmap)) {
    geocode_OSM("India")
    geocode_OSM("CBS Weg 1, Heerlen")
    geocode_OSM("CBS Weg 1, Heerlen", projection = 28992)

    data(metro)

    # sample 5 cities from the metro dataset
    five_cities <- metro[sample(length(metro), 5), ]

    # obtain geocode locations from their long names
    five_cities_geocode <- geocode_OSM(five_cities$name_long, as.sf = TRUE)

    # change to interactive mode
    current.mode <- tmap_mode("view")

    # plot metro coordinates in red and geocode coordinates in blue
    # zoom in to see the differences
    tm_shape(five_cities) +
    	tm_dots(col = "blue") +
    tm_shape(five_cities_geocode) +
    	tm_dots(col = "red")

    # restore current mode
    tmap_mode(current.mode)
}

## End(Not run)

Get aspect ratio

Description

Get the aspect ratio of a shape object, a tmap object, or a bounding box

Usage

get_asp_ratio(x, is.projected = NA, width = 700, height = 700, res = 100)

Arguments

Details

The arguments width, height, and res are passed on to png. If x is a tmap object, a temporarily png image is created to calculate the aspect ratio of a tmap object. The default size of this image is 700 by 700 pixels at 100 dpi.

Value

aspect ratio

Examples

if (require(tmap) && packageVersion("tmap") >= "2.0") {
    data(World)

    get_asp_ratio(World)

    get_asp_ratio(bb(World))

    tm <- qtm(World)
    get_asp_ratio(tm)
}

## Not run: 
    get_asp_ratio("Germany") #note: bb("Germany") uses geocode_OSM("Germany")

## End(Not run)

Get neighbours list from spatial objects

Description

Get neighbours list from spatial objects. The output is similar to the function poly2nb of the spdep package, but uses sf instead of sp.

Usage

get_neighbours(x)

Arguments

Value

A list where the items correspond to the features. Each item is a vector of neighbours.

Map coloring

Description

Color the polygons of a map such that adjacent polygons have different colors. This function returns the color indices

Usage

map_coloring(x, algorithm = "greedy", ncols = 8, minimize = FALSE, ...)

Arguments

Details

As of tmaptools 3.3, the deprecated color functions get_brewer_pal and palette_explorer, have been removed. These have been replaced c4a and c4a_gui respectively from the package cols4all. Therefore, map_coloring will return color indices and will ignore the input arguments palette and contrast. See example.

Value

A vector of color indices.

Examples

if (require(tmap) && require(cols4all)) {
    data(World)

    ## using cols4all directly
    indices <- map_coloring(World)
    pal <- c4a("brewer.set2", n = max(indices))
    World$color = pal[indices]
    tm_shape(World) +
        tm_polygons("color", fill.scale = tm_scale_asis()) +
        tm_crs("auto")

    # using map_coloring via "MAP_COLORS" in tmap
    tm_shape(World) +
        tm_polygons("MAP_COLORS", tm_scale(values = "brewer.set2")) +
        tm_crs("auto")

    # other example
    data(NLD_prov, NLD_muni)
    tm_shape(NLD_prov) +
    	tm_fill("name",
    	        fill.legend = tm_legend_hide()) +
    tm_shape(NLD_muni) +
    	tm_polygons("MAP_COLORS",
    	            fill_alpha = .25,
    	            fill.scale = tm_scale(values = "brewer.greys")) +
    tm_shape(NLD_prov) +
    	tm_borders(lwd=2) +
    	tm_text("name", options = opt_tm_text(shadow = TRUE)) +
    tm_title("Dutch provinces and\nmunicipalities", bg.color="white")
}

Read GPX file

Description

Read a GPX file. By default, it reads all possible GPX layers, and only returns shapes for layers that have any features.

Usage

read_GPX(
  file,
  layers = c("waypoints", "routes", "tracks", "route_points", "track_points"),
  remove.empty.layers = TRUE,
  as.sf = TRUE
)

Arguments

Details

Note that this function returns sf objects, but still uses methods from sp and rgdal internally.

Value

a list of sf objects, one for each layer

Read Open Street Map data

Description

Read Open Street Map data. OSM tiles are read and returned as a spatial raster. Vectorized OSM data is not supported anymore (see details).

Usage

read_osm(
  x,
  zoom = NULL,
  type = "osm",
  minNumTiles = NULL,
  mergeTiles = NULL,
  use.colortable = FALSE,
  ...
)

Arguments

Details

As of version 2.0, read_osm cannot be used to read vectorized OSM data anymore. The reason is that the package that was used under the hood, osmar, has some limitations and is not actively maintained anymore. Therefore, we recommend the package osmdata. Since this package is very user-friendly, there was no reason to use read_osm as a wrapper for reading vectorized OSM data.

Value

The output of read_osm is a raster object.

Examples

## Not run: 
if (require(tmap)) {
    #### Choropleth with OSM background

    # load Netherlands shape
    data(NLD_muni)

    # read OSM raster data
    osm_NLD <- read_osm(NLD_muni, ext=1.1)

    # plot with regular tmap functions
    tm_shape(osm_NLD) +
    	tm_rgb() +
    tm_shape(NLD_muni) +
    	tm_polygons("population", convert2density=TRUE, style="kmeans", alpha=.7, palette="Purples")

    #### A close look at the building of Statistics Netherlands in Heerlen

    # create a bounding box around the CBS (Statistics Netherlands) building
    CBS_bb <- bb("CBS Weg 11, Heerlen", width=.003, height=.002)

    # read Microsoft Bing satellite and OpenCycleMap OSM layers
    CBS_osm1 <- read_osm(CBS_bb, type="bing")
    CBS_osm2 <- read_osm(CBS_bb, type="opencyclemap")

    # plot OSM raster data
    qtm(CBS_osm1)
    qtm(CBS_osm2)

}

## End(Not run)

Reverse geocodes a location using OpenStreetMap Nominatim

Description

Reverse geocodes a location (based on spatial coordinates) to an address. It uses OpenStreetMap Nominatim. For processing large amount of queries, please read the usage policy (https://operations.osmfoundation.org/policies/nominatim/).

Usage

rev_geocode_OSM(
  x,
  y = NULL,
  zoom = NULL,
  projection = 4326,
  as.data.frame = NA,
  server = "https://nominatim.openstreetmap.org",
  params = NULL
)

Arguments

Value

A data frame or a list with all attributes that are contained in the search result

See Also

geocode_OSM

Examples

## Not run: 
if (require(tmap)) {
    data(metro)

    # sample five cities from metro dataset
    set.seed(1234)
    five_cities <- metro[sample(length(metro), 5), ]

    # obtain reverse geocode address information
    addresses <- rev_geocode_OSM(five_cities, zoom = 6)
    five_cities <- sf::st_sf(data.frame(five_cities, addresses))

    # change to interactive mode
    current.mode <- tmap_mode("view")
    tm_shape(five_cities) +
    	tm_markers(text="name")

    # restore current mode
    tmap_mode(current.mode)
}

## End(Not run)

Simplify shape

Description

Simplify a shape consisting of polygons or lines. This can be useful for shapes that are too detailed for visualization, especially along natural borders such as coastlines and rivers. The number of coordinates is reduced.

Usage

simplify_shape(shp, fact = 0.1, keep.units = FALSE, keep.subunits = FALSE, ...)

Arguments

Details

This function is a wrapper of ms_simplify. In addition, the data is preserved. Also sf objects are supported.

Value

sf object

Examples

## Not run: 
if (require(tmap)) {
    data(World)

    # show different simplification factors
    tm1 <- qtm(simplify_shape(World, fact = 0.05), title="Simplify 0.05")
    tm2 <- qtm(simplify_shape(World, fact = 0.1), title="Simplify 0.1")
    tm3 <- qtm(simplify_shape(World, fact = 0.2), title="Simplify 0.2")
    tm4 <- qtm(simplify_shape(World, fact = 0.5), title="Simplify 0.5")
    tmap_arrange(tm1, tm2, tm3, tm4)

    # show different options for keeping smaller (sub)units
    tm5 <- qtm(simplify_shape(World, keep.units = TRUE, keep.subunits = TRUE),
        title="Keep units and subunits")
    tm6 <- qtm(simplify_shape(World, keep.units = TRUE, keep.subunits = FALSE),
        title="Keep units, ignore small subunits")
    tm7 <- qtm(simplify_shape(World, keep.units = FALSE),
        title="Ignore small units and subunits")
    tmap_arrange(tm5, tm6, tm7)
}

## End(Not run)