%\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{Overview of mapping with mapmisc} \documentclass[12pt]{article} \usepackage[margin=1in]{geometry} \usepackage{graphicx,hyperref} \DeclareGraphicsExtensions{.png} \usepackage{caption,subcaption} \providecommand{\subfloat}[2][need a sub-caption]{\subcaptionbox{#1}{#2}} \title{The mapmisc package} \author{Patrick Brown} \begin{document} \maketitle <>= library('knitr') opts_chunk$set(echo=TRUE,message=FALSE, dev='png', fig.align='center',fig.pos='ht',width=80, fig.height=5,fig.width=3,out.width="0.3\\textwidth", dev.args=list(bg='white') ) options("rgdal_show_exportToProj4_warnings"="none") library('mapmisc') haveRgdal = requireNamespace('rgdal', quietly=TRUE) if(!exists('mapmiscCachePath')) mapmiscCachePath = system.file('extdata', package='mapmisc') if(!exists('mapmiscCacheReadOnly')) mapmiscCacheReadOnly = TRUE if(!exists('fact')) fact = 0.6 # uncomment for nicer images # fact = 1 options( mapmiscCachePath = mapmiscCachePath, mapmiscCacheReadOnly = mapmiscCacheReadOnly ) opts_chunk$set(dpi=72*fact) @ <>= if(!exists('testing')) testing=FALSE if(testing) { cat('cache:', getOption('mapmiscCachePath')) cat('\nread only:', getOption('mapmiscCacheReadOnly')) cat('\n') options(mapmiscVerbose = TRUE) } @ This document will be incomplete if {\tt rgdal} is unavailable or there is on internet connection when this document is compiled. The full document is at \url{http://diseasemapping.r-forge.r-project.org/}. \section{Introduction} This package provides a few utilities for making nice maps in R, with an emphasis on enabling maps in short tidy code chunks which are suitable for Sweave and knitr documents. The package duplicates the capabilities of packages such as classInt, geonames, and OpenStreetMap, and the price of having tidier code is much of the flexibility from these other packages has been lost here. The meuse data from the sp package <>= library('sp') data('meuse') coordinates(meuse) <- ~x+y class(meuse) @ <>= library('rgdal') proj4string(meuse) <- CRS("+init=epsg:28992") meuseLL = spTransform(meuse, CRS("+init=epsg:4326")) @ <>= if(haveRgdal) { proj4string(meuse) <- CRS("+init=epsg:28992") meuseLL = spTransform(meuse, CRS("+init=epsg:4326")) } else { proj4string(meuse) <- CRS( "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs" ) meuseLL = new("SpatialPointsDataFrame" , data = structure(list(cadmium = c(11.7, 8.6, 6.5, 2.6, 2.8, 3, 3.2, 2.8, 2.4, 1.6, 1.4, 1.8, 11.2, 2.5, 2, 9.5, 7, 7.1, 8.7, 12.9, 5.5, 2.8, 2.9, 1.7, 1.4, 1.5, 1.3, 1.3, 1.8, 1.5, 1.3, 1.3, 2.1, 1.2, 1.6, 3, 5.8, 7.9, 8.1, 14.1, 2.4, 1.2, 1.7, 1.3, 4.2, 4.3, 3.1, 1.7, 1.7, 2.1, 1.7, 8.2, 17, 12, 9.4, 8.2, 2.6, 3.5, 10.9, 7.3, 9.4, 8.3, 7, 6.8, 7.4, 6.6, 7.8, 0.7, 3.4, 3.9, 3.5, 4.7, 3.9, 3.1, 2.9, 3.9, 2.7, 2.5, 5.6, 9.4, 10.8, 18.1, 6.3, 2.1, 1.8, 2.2, 2.8, 2.4, 2.6, 1.8, 2, 2.7, 2.7, 0.9, 0.9, 0.4, 0.8, 0.4, 0.4, 0.4, 0.8, 0.4, 0.4, 0.4, 0.2, 0.2, 0.4, 0.2, 0.2, 0.4, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 3.7, 0.2, 0.2, 0.2, 0.2, 1.2, 2.4, 0.2, 0.2, 0.2, 0.2, 0.4, 1.6, 0.8, 1.2, 0.2, 0.2, 0.8, 1, 0.8, 1.2, 2, 1.5, 1.1, 0.8, 0.8, 2.1, 2.5, 3.8, 3.2, 3.1, 1.5, 1.2, 0.8, 3.1, 2.1, 0.8, 2.7), copper = c(85, 81, 68, 81, 48, 61, 31, 29, 37, 24, 25, 25, 93, 31, 27, 86, 74, 69, 69, 95, 53, 35, 35, 24, 26, 22, 27, 24, 22, 21, 29, 21, 27, 26, 27, 32, 50, 67, 77, 108, 32, 21, 22, 21, 51, 50, 38, 26, 24, 32, 24, 47, 128, 117, 104, 76, 36, 34, 90, 80, 78, 77, 65, 66, 72, 75, 75, 22, 55, 47, 46, 55, 49, 39, 45, 48, 36, 36, 68, 88, 85, 76, 63, 32, 25, 27, 36, 41, 33, 27, 38, 37, 34, 19, 22, 26, 25, 22, 24, 18, 23, 22, 23, 23, 23, 18, 20, 22, 27, 25, 30, 23, 24, 21, 24, 21, 23, 53, 33, 22, 23, 25, 30, 47, 31, 20, 16, 17, 26, 34, 18, 23, 20, 14, 46, 29, 19, 31, 27, 29, 22, 20, 20, 29, 38, 39, 35, 42, 24, 20, 20, 39, 31, 21, 27), lead = c(299, 277, 199, 116, 117, 137, 132, 150, 133, 80, 86, 97, 285, 183, 130, 240, 133, 148, 207, 284, 194, 123, 110, 85, 75, 76, 73, 67, 87, 65, 78, 64, 77, 80, 82, 97, 166, 217, 219, 405, 102, 48, 65, 62, 281, 294, 211, 135, 112, 162, 94, 191, 405, 654, 482, 276, 180, 207, 541, 310, 210, 158, 141, 144, 181, 173, 399, 45, 325, 268, 252, 315, 260, 237, 228, 241, 201, 204, 429, 462, 333, 464, 159, 116, 81, 131, 216, 145, 163, 129, 148, 214, 226, 54, 70, 73, 87, 76, 81, 68, 66, 49, 63, 48, 51, 50, 39, 48, 64, 84, 67, 49, 48, 56, 65, 84, 75, 250, 81, 72, 86, 94, 244, 297, 96, 56, 49, 50, 148, 162, 37, 80, 80, 49, 42, 48, 41, 73, 146, 95, 72, 51, 54, 136, 170, 179, 200, 258, 93, 68, 49, 173, 119, 51, 124), zinc = c(1022, 1141, 640, 257, 269, 281, 346, 406, 347, 183, 189, 251, 1096, 504, 326, 1032, 606, 711, 735, 1052, 673, 402, 343, 218, 200, 194, 207, 180, 240, 180, 208, 198, 250, 192, 213, 321, 569, 833, 906, 1454, 298, 167, 176, 258, 746, 746, 464, 365, 282, 375, 222, 812, 1548, 1839, 1528, 933, 432, 550, 1571, 1190, 907, 761, 659, 643, 801, 784, 1060, 119, 778, 703, 676, 793, 685, 593, 549, 680, 539, 560, 1136, 1383, 1161, 1672, 765, 279, 241, 317, 545, 505, 420, 332, 400, 553, 577, 155, 224, 180, 226, 186, 198, 187, 199, 157, 203, 143, 136, 117, 113, 130, 192, 240, 221, 140, 128, 166, 191, 232, 203, 722, 210, 198, 139, 253, 703, 832, 262, 142, 119, 152, 415, 474, 126, 210, 220, 133, 141, 158, 129, 206, 451, 296, 189, 154, 169, 403, 471, 612, 601, 783, 258, 214, 166, 496, 342, 162, 375), elev = c(7.909, 6.983, 7.8, 7.655, 7.48, 7.791, 8.217, 8.49, 8.668, 9.049, 9.015, 9.073, 7.32, 8.815, 8.937, 7.702, 7.16, 7.1, 7.02, 6.86, 8.908, 8.99, 8.83, 9.02, 8.976, 8.973, 8.507, 8.743, 8.973, 9.043, 8.688, 8.727, 8.328, 7.971, 7.809, 6.986, 7.756, 7.784, 7, 6.92, 7.516, 8.18, 8.694, 9.28, 7.94, 6.36, 7.78, 8.18, 9.42, 8.867, 8.292, 8.06, 7.98, 7.9, 7.74, 7.552, 7.76, 6.74, 6.68, 7.4, 7.44, 7.36, 7.2, 7.22, 7.36, 5.18, 5.8, 7.64, 6.32, 5.76, 6.48, 6.48, 6.32, 6.32, 6.16, 6.56, 6.9, 7.54, 6.42, 6.28, 6.34, 7.307, 5.7, 7.72, 7.932, 7.82, 8.575, 8.536, 8.504, 8.659, 7.633, 8.538, 7.68, 7.564, 7.76, 7.653, 7.951, 8.176, 8.468, 8.41, 7.61, 7.792, 8.76, 9.879, 9.097, 9.095, 9.717, 9.924, 9.404, 10.52, 8.84, 8.472, 9.634, 9.206, 8.47, 8.463, 9.691, 8.704, 9.42, 9.573, 9.555, 8.779, 8.54, 8.809, 9.523, 9.811, 9.604, 9.732, 9.518, 9.72, 9.036, 9.528, 9.155, 10.08, 9.97, 10.136, 10.32, 9.041, 7.86, 8.741, 7.822, 7.78, 8.121, 8.231, 8.351, 7.3, 7.536, 7.706, 8.07, 8.226, 8.128, 8.577, 8.429, 9.406, 8.261), dist = c(0.00135803, 0.0122243, 0.103029, 0.190094, 0.27709, 0.364067, 0.190094, 0.0921516, 0.184614, 0.309702, 0.315116, 0.228123, 0, 0.113932, 0.168336, 0, 0.0122243, 0.0122243, 0, 0, 0.0703468, 0.0975136, 0.113932, 0.342321, 0.385804, 0.429289, 0.315116, 0.320574, 0.249863, 0.489064, 0.472778, 0.287957, 0.271622, 0.385807, 0.418417, 0.244474, 0.135709, 0.0484965, 0, 0.00135803, 0.135709, 0.26622, 0.211843, 0.320572, 0.081222, 0.190086, 0.287941, 0.423826, 0.554289, 0.603225, 0.614071, 0.00135803, 0, 0.0054321, 0.0054321, 0.0054321, 0.146578, 0.135684, 0.0703333, 0.0484831, 0.0054321, 0.0054321, 0.0316663, 0.0122243, 0.0122243, 0.0373395, 0.211846, 0.451037, 0.575877, 0.0756869, 0.12481, 0.103024, 0.157469, 0.200976, 0.200976, 0.26622, 0.298835, 0.0812247, 0.070355, 0.0122243, 0, 0.0537723, 0.0593662, 0.211843, 0.12481, 0.12481, 0.0921516, 0.113941, 0.179216, 0.233596, 0.336861, 0.070355, 0.0593662, 0.255341, 0.364067, 0.429295, 0.380328, 0.249874, 0.266212, 0.451037, 0.30971, 0.293359, 0.532351, 0.619513, 0.684725, 0.809742, 0.880389, 0.749591, 0.575752, 0.581484, 0.49452, 0.32058, 0.336851, 0.249852, 0.0756869, 0.0756869, 0.162853, 0.0974916, 0.440142, 0.4619, 0.222701, 0.103024, 0.0921353, 0.0484884, 0.168331, 0.38581, 0.489064, 0.57602, 0.0812194, 0.0373369, 0.771698, 0.336829, 0.228123, 0.597761, 0.44558, 0.396675, 0.581478, 0.287966, 0.233596, 0.364067, 0.331454, 0.211846, 0.103029, 0.070355, 0.146576, 0.0537723, 0.119286, 0.0921435, 0.249863, 0.37494, 0.423837, 0.423837, 0.27709, 0.358606, 0.0122243), om = c(13.6, 14, 13, 8, 8.7, 7.8, 9.2, 9.5, 10.6, 6.3, 6.4, 9, 15.4, 8.4, 9.1, 16.2, 16, 16, 13.7, 14.8, 10.2, 7.2, 7.2, 7, 6.9, 6.3, 5.6, 4.4, 5.3, 4.8, 2.6, 1, 2.4, 1.9, 3.1, 1.6, 3.5, 8.1, 7.9, 9.5, 1.4, NA, NA, 2, 5.1, 5.3, 4.5, 4.9, 4.5, 5.5, 3.4, 11.1, 12.3, 16.5, 13.9, 8.1, 3.1, 5.8, 10.2, 12, 14.1, 14.5, 14.8, 13.3, 15.2, 11.4, 9, 3.6, 6.9, 7, 6.2, 6.5, 5.7, 7, 7.3, 8.2, 4.3, 4.4, 8.2, 8.5, 9.6, 17, 12.8, 5.9, 2.9, 4.5, 10.7, 9.4, 9, 7, 6.5, 9.4, 10.2, 6.4, 7.6, 7, 5.6, 6.5, 6.6, 5.9, 6.5, 6.4, 7.2, 6.6, 4.3, 5.3, 4.1, 6.1, 7.5, 8.8, 5.7, 6.1, 7.1, 4.1, 6, 6.6, 6.8, 9.1, 5.9, 4.9, 7.1, 8.1, 8.3, 10, 5.9, 5, 4.5, 5.4, 2.3, 7.5, 4.6, 5.8, 4.4, 4.4, 4.5, 5.2, 4.6, 6.9, 7, 5.4, 6.2, 5, 5.1, 8.1, 8, 8.8, 9.3, 8.4, 7.7, 5.7, 4.7, 9.1, 6.5, 5.7, 5.5), ffreq = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor"), soil = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 3L, 2L, 1L, 2L, 2L, 2L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L), .Label = c("1", "2", "3"), class = "factor"), lime = structure(c(2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("0", "1"), class = "factor"), landuse = structure(c(4L, 4L, 4L, 11L, 4L, 11L, 4L, 2L, 2L, 15L, 9L, 3L, 15L, 4L, 4L, 15L, 15L, 15L, 15L, NA, 5L, 5L, 3L, 4L, 15L, 15L, 2L, 3L, 4L, 3L, 6L, 3L, 4L, 6L, 6L, 2L, 2L, 5L, 15L, 15L, 5L, 11L, 15L, 4L, 4L, 5L, 4L, 4L, 7L, 7L, 2L, 4L, 15L, 15L, 15L, 15L, 10L, 10L, 10L, 15L, 15L, 10L, 15L, 4L, 15L, 15L, 4L, 10L, 7L, 2L, 2L, 15L, 15L, 4L, 15L, 15L, 4L, 5L, 15L, 15L, 15L, 15L, 15L, 15L, 4L, 15L, 15L, 15L, 4L, 5L, 5L, 5L, 5L, 15L, 15L, 5L, 5L, 5L, 4L, 15L, 15L, 12L, 15L, 5L, 4L, 10L, 4L, 4L, 10L, 13L, 8L, 10L, 4L, 4L, 4L, 15L, 13L, 7L, 4L, 1L, 15L, 14L, 10L, 4L, 4L, 4L, 5L, 5L, 5L, 15L, 4L, 15L, 1L, 5L, 5L, 5L, 15L, 15L, 15L, 4L, 4L, 10L, 15L, 7L, 7L, 15L, 15L, 4L, 5L, 4L, 5L, 4L, 4L, 15L, 15L), .Label = c("Aa", "Ab", "Ag", "Ah", "Am", "B", "Bw", "DEN", "Fh", "Fw", "Ga", "SPO", "STA", "Tv", "W"), class = "factor"), dist.m = c(50, 30, 150, 270, 380, 470, 240, 120, 240, 420, 400, 300, 20, 130, 220, 10, 10, 10, 10, 10, 80, 140, 160, 440, 490, 530, 400, 400, 330, 630, 570, 390, 360, 500, 550, 340, 210, 60, 10, 20, 170, 320, 260, 360, 100, 200, 320, 480, 660, 690, 710, 10, 10, 10, 10, 20, 200, 140, 70, 20, 10, 10, 20, 10, 20, 20, 270, 560, 750, 80, 130, 110, 200, 260, 270, 320, 360, 80, 100, 70, 20, 50, 80, 290, 160, 160, 140, 150, 220, 280, 450, 70, 70, 340, 470, 630, 460, 320, 320, 540, 420, 380, 500, 760, 750, 1000, 860, 680, 500, 650, 630, 410, 390, 310, 70, 70, 150, 80, 450, 530, 240, 70, 70, 60, 190, 450, 550, 650, 100, 170, 860, 410, 290, 680, 540, 480, 720, 380, 310, 430, 370, 290, 150, 70, 220, 80, 120, 120, 260, 440, 540, 520, 350, 460, 40)), .Names = c("cadmium", "copper", "lead", "zinc", "elev", "dist", "om", "ffreq", "soil", "lime", "landuse", "dist.m"), row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "69", "75", "76", "79", "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "123", "160", "163", "70", "71", "91", "92", "93", "94", "95", "96", "97", "98", "99", "100", "101", "102", "103", "104", "105", "106", "108", "109", "110", "111", "112", "113", "114", "115", "116", "117", "118", "119", "120", "121", "122", "124", "125", "126", "127", "128", "129", "130", "131", "132", "133", "134", "135", "136", "161", "162", "137", "138", "140", "141", "142", "143", "144", "145", "146", "147", "148", "149", "150", "151", "152", "153", "154", "155", "156", "157", "158", "159", "164"), class = "data.frame") , coords.nrs = numeric(0) , coords = structure(c(5.75853624181867, 5.75786301834362, 5.75985541797357, 5.76174577016607, 5.76186270784184, 5.76303964144362, 5.75984329328111, 5.75787741999485, 5.75833785418251, 5.76078279124001, 5.76019504951489, 5.75792464994252, 5.75569670452327, 5.75704339874017, 5.75763497327733, 5.75506340377509, 5.75409907584572, 5.75310703264216, 5.75211553795547, 5.75110878811622, 5.75234773689981, 5.75323125695777, 5.75322747806444, 5.75964010271363, 5.75954726088018, 5.75982879999144, 5.75756454334514, 5.75705965370699, 5.75625272368661, 5.7624755086065, 5.75932955716756, 5.75569265030061, 5.75499212770497, 5.75676836289353, 5.75679101010027, 5.75328909953012, 5.75218663874741, 5.75074082096063, 5.75000319333192, 5.74864324807759, 5.75021349112691, 5.75115778755962, 5.74959027645615, 5.74899635068879, 5.74823209944538, 5.74732509154481, 5.74718691918422, 5.74716195489008, 5.74698031276162, 5.74596136608848, 5.74504748711013, 5.74655555462392, 5.74464376906767, 5.7427898090949, 5.74069409765227, 5.73871109250081, 5.7401974707017, 5.74176905188958, 5.74459520513393, 5.73752233818859, 5.73640339232072, 5.73595074325204, 5.7358560163478, 5.73477883018096, 5.73363126295385, 5.73249973693215, 5.7355510507878, 5.73869821749255, 5.74096673840647, 5.73084191790459, 5.73066521079125, 5.72976821070395, 5.72893346997607, 5.73040169842579, 5.7292838718135, 5.73018003212326, 5.7296449629613, 5.72612524059285, 5.72758463522269, 5.72857646604616, 5.72993420415961, 5.74519169091642, 5.73186737889372, 5.75916262406516, 5.73562013746168, 5.73683830959438, 5.72915670865268, 5.72816580082885, 5.72704054976721, 5.72600468925749, 5.7292244546065, 5.72319018528861, 5.72456653418406, 5.73658732680482, 5.73308151854868, 5.73458799084597, 5.73541282448076, 5.73300925888455, 5.73139490487106, 5.73174407181168, 5.73761241300685, 5.73995238255195, 5.74278684307447, 5.74405388185952, 5.74542037012928, 5.74954073570442, 5.74776752152038, 5.74701360362797, 5.74469637173753, 5.74138719979516, 5.73986467036325, 5.73821482063726, 5.73927477424651, 5.73908890124233, 5.73522952314037, 5.73634968875584, 5.73806938376514, 5.74761229151849, 5.74545176360193, 5.74355004601903, 5.74026068485757, 5.73805628328566, 5.74101641782673, 5.74695046270072, 5.74468707713721, 5.74973080115803, 5.74994977613158, 5.74807608601104, 5.755509802491, 5.75318606019536, 5.74595972853237, 5.74561437991626, 5.75636110555308, 5.74978412489021, 5.74192818602271, 5.74057003898231, 5.742963927314, 5.73053201209076, 5.7293067454292, 5.73002413981022, 5.73214744173562, 5.73528038404125, 5.73660928002864, 5.73510587704988, 5.73357307132981, 5.73225323517617, 5.72910928823648, 5.72572798800468, 5.7306995133442, 5.72921791366753, 5.73141645736525, 5.73001472930627, 5.72702698674827, 5.73544365891206, 5.7519552418591, 50.99156215511, 50.9910879013154, 50.9908927417686, 50.990410234025, 50.9890255693492, 50.9883925350409, 50.989391639099, 50.9893350252852, 50.9881470197354, 50.9875728602917, 50.9870983443726, 50.9863505688173, 50.9891262510716, 50.9883399268111, 50.9875200489596, 50.9882922993244, 50.9870189386853, 50.9858355464053, 50.9847150661029, 50.9834597923402, 50.9834648814659, 50.9835429922083, 50.9830665940866, 50.985392244576, 50.9844756781641, 50.9840702734161, 50.9840685476876, 50.9832611628361, 50.9838570025257, 50.9855718580696, 50.9822021933378, 50.9814857242104, 50.9811373851867, 50.9806732980246, 50.9799451272401, 50.9800641170388, 50.9811013168348, 50.9819418293055, 50.9823036972055, 50.9813910851934, 50.980073798991, 50.978839360128, 50.9786824764636, 50.9773899367106, 50.9798013382546, 50.9785277433723, 50.9772427651041, 50.9758675461797, 50.9744838253139, 50.973444278356, 50.9730965404907, 50.9802739739873, 50.9797944919884, 50.9794226673543, 50.9790875067904, 50.9785901667863, 50.9774530485726, 50.9781493785474, 50.9790395771904, 50.9776769051776, 50.9765746548795, 50.975083868009, 50.9738526792524, 50.9726064660439, 50.9714593310135, 50.9705458472601, 50.9692872556725, 50.9675518827339, 50.9663584436644, 50.9677732219455, 50.9669827267056, 50.9669584176528, 50.9657294109239, 50.9660037225325, 50.9649912897607, 50.964907738672, 50.9637767264194, 50.9651893455328, 50.966200804956, 50.967501269146, 50.9682253554512, 50.9623184459304, 50.9697477228722, 50.987991563045, 50.9708960539211, 50.9719620563834, 50.9670411250794, 50.9658665071403, 50.9638383338299, 50.9623402307929, 50.9627352619485, 50.9628607980544, 50.9642141127235, 50.961373935548, 50.9607911980511, 50.9642743650748, 50.9660606704467, 50.9663465290704, 50.9655693088616, 50.9627637462787, 50.9632854648137, 50.9640154455151, 50.9661011798284, 50.9698006936087, 50.9716122247639, 50.9723724467692, 50.969546489231, 50.9677960057513, 50.9663649654033, 50.9673549219818, 50.9674404693051, 50.9695039323028, 50.9703546560202, 50.9721170381136, 50.971957914276, 50.9732399491374, 50.9747179131779, 50.9624277909682, 50.975657117018, 50.9746112877967, 50.9764101477016, 50.9767404469175, 50.9784123576463, 50.9797873490838, 50.9780325406988, 50.976632574868, 50.9755082807005, 50.9744174934376, 50.9871223274747, 50.9850262984647, 50.9695700824649, 50.9764296562356, 50.9849443104351, 50.974358224509, 50.9710656664222, 50.9703237436441, 50.9687613395114, 50.9643853416156, 50.9642361580936, 50.9630834621507, 50.9595984667859, 50.9598227996314, 50.9586322668089, 50.957441276325, 50.9580121708692, 50.9566138504223, 50.9567939672551, 50.9576039129667, 50.9582005102161, 50.9599307576448, 50.9608231273358, 50.9618160610117, 50.9619956579452, 50.9625998664103, 50.9608321490479), .Dim = c(155L, 2L), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "69", "75", "76", "79", "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "123", "160", "163", "70", "71", "91", "92", "93", "94", "95", "96", "97", "98", "99", "100", "101", "102", "103", "104", "105", "106", "108", "109", "110", "111", "112", "113", "114", "115", "116", "117", "118", "119", "120", "121", "122", "124", "125", "126", "127", "128", "129", "130", "131", "132", "133", "134", "135", "136", "161", "162", "137", "138", "140", "141", "142", "143", "144", "145", "146", "147", "148", "149", "150", "151", "152", "153", "154", "155", "156", "157", "158", "159", "164"), c("x", "y"))) , bbox = structure(c(5.72319018528861, 50.9566138504223, 5.76303964144362, 50.99156215511), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y" ), c("min", "max"))) , proj4string = new("CRS" , projargs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" ) ) } @ Get elevation data with <>= library('raster') nldElev = raster::getData("alt", lon=raster::xmin(meuseLL), lat=raster::ymin(meuseLL), path=tempdir() ) raster::getData("SRTM", lon=raster::xmin(meuseLL), lat=raster::ymin(meuseLL), path=tempdir() ) nldElev = raster(file.path(tempdir(), "srtm_38_02.tif")) nldElev = crop(nldElev, extent(meuseLL)) nldElev = projectRaster(nldElev, crs=proj4string(meuse)) @ Or simply do <>= library('mapmisc') data('netherlands') @ The elevation data is a Raster. <>= class(nldElev) nldElev = crop(nldElev, extend(extent(meuse), 1000)) @ \section{Downloading background maps and city locations} Get a background map covering the extent of the meuse data <>= library('mapmisc') if(haveRgdal) nldTiles = openmap(meuse) @ {\tt nldTiles} is a Raster with the same projection as {\tt meuse} <>= class(nldTiles) projection(nldTiles) projection(meuse) @ Maps which can be downloaded are shown at \url{http://diseasemapping.r-forge.r-project.org/openmap/} {\tt GNcities} is a wrapper for the function of the same name in the geonames package. <>= nldCities = GNcities(meuse, maxRows=6) @ A SpatialPointsDataFrame, with same map projection. <>= class(nldCities) names(nldCities) projection(nldCities) @ \section{Making maps} The {\tt map.new} function sets up a map in the current plot window with the correct limits and aspect ratio for the object supplied, and without margins or white space. {\tt scaleBar} adds a scale and north arrow. It uses the map projection of the argument supplied to calculate distances and find north. Some of the code in the following sections require the {\tt rgdal} package to run. <>= # plot the data locations map.new(meuse) plot(nldTiles, add=TRUE) points(meuse,col="red", cex=0.3) scaleBar(crs=proj4string(meuse),pos="topleft", bg="white") # plot city names map.new(meuse) plot(nldTiles, add=TRUE) points(nldCities) text(nldCities, labels=nldCities$name, pos=3) scaleBar(proj4string(meuse),pos="topleft", bg="white") # plot elevation map.new(meuse, legendRight=TRUE) plot(nldTiles, add=TRUE) plot(nldElev,add=TRUE,col=terrain.colors(8),alpha=0.6,legend.mar=2, legend.line=0,) scaleBar(meuse,pos="topleft",bg="white") @ \section{Legends} Create a colour scale for plotting copper concentrations <>= cuScale = colourScale(meuse$copper, breaks=5, style='equal', opacity=0.8, dec=-1, firstBreak=0) @ and elevation, with transparency decreasing as elevation increases. <>= elevScale = colourScale(nldElev, style='equal', breaks=6, col=terrain.colors, firstBreak=0, dec=-1,opacity=c(0.2, 0.9)) @ Soil type is a categorical variable, create a factor and create a colour scale of unique values <>= meuse$soilFac = factor(meuse$soil, levels=c(1,2,3), labels=c("Calcareous","Non-Calc's","Red Brick")) soilScale = colourScale(meuse$soilFac, col="Set2") @ <>= map.new(meuse) plot(nldTiles, add=TRUE) plot(meuse, col=cuScale$plot,add=TRUE,pch=16) legendBreaks("bottomright", breaks=cuScale, title="gals/firkin") map.new(meuse) plot(nldTiles, add=TRUE) plot(meuse, col=soilScale$plot,add=TRUE,pch=16) legendBreaks("bottomright", breaks=soilScale, title="soil type", cex=0.7,bg="white") map.new(meuse) plot(nldTiles, add=TRUE) image(nldElev, breaks=elevScale$breaks, col=elevScale$colOpacity, legend=FALSE,add=TRUE) legendBreaks("left", breaks=elevScale, title='Metres',bg="white") @ \section{More plots} Rotate the data 50 degrees clockwise with an oblique mercator projection. <>= if(requireNamespace('rgdal', quietly=TRUE)){ meuseRot = spTransform(meuse, omerc(meuse, -50)) tilesRot = openmap(meuseRot, fact=2) elevRot = projectRaster(nldElev, crs=projection(meuseRot)) nldCitiesRot = spTransform(nldCities, CRS(projection(meuseRot))) } @ <>= if(!requireNamespace('rgdal', quietly=TRUE)) { tilesRot = nldTiles elevRot = nldElev meuseRot = meuse nldCitiesRot = nldCities } @ And create new plots <>= # first elevation map.new(meuseRot) plot(tilesRot, add=TRUE) plot(elevRot,add=TRUE,alpha=0.5,col=terrain.colors(8), legend=FALSE) points(nldCitiesRot) text(nldCitiesRot, labels=nldCitiesRot$name, pos=3) scaleBar(meuseRot,pos="topleft", bg="white") # then data locations map.new(meuseRot) plot(tilesRot, add=TRUE) points(meuseRot,col="red", cex=0.3) scaleBar(proj4string(meuseRot), bg="white") @ \section{Inset maps} <>= world = openmap( extent(-10,30,40,60), path="osm-no-labels") @ <>= # not rotated map.new(meuse,legendRight=TRUE) plot(nldTiles, add=TRUE) points(meuse) if(requireNamespace('rgdal', quietly=TRUE)) { scaleBar(proj4string(meuse),pos="bottomright", bg="white") insetMap(crs=meuse, pos="topright",map=world) } # rotated map.new(meuseRot) plot(tilesRot, add=TRUE) points(meuseRot,col="red", cex=0.3) if(requireNamespace('rgdal', quietly=TRUE)) { scaleBar(proj4string(meuseRot), bg="white") insetMap(meuseRot, "bottomleft",map=world) } @ \end{document}