• Artículos
  • Documentos
    • Libro de Energía Solar Fotovoltaica
  • Omnia sunt Communia!

Omnia sunt Communia!

~ Sobre software, documentación y ciencia libres

Omnia sunt Communia!

Archivos de etiqueta: lattice

Mapping Flows in R … with data.table and lattice

14 martes Abr 2015

Posted by Oscar Perpiñán Lamigueiro in R-english, visualization

≈ 6 comentarios

Etiquetas

data.table, lattice, R, spatial data

Some days ago James Cheshire published the post Mapping Flows in R. I have implemented an alternative (faster) version using data.table to read and join the datasets (and lattice to display the results). If you are new to data.table you should read this wiki and this cheatsheet.

This is the code:


### DATA SECTION
library(data.table)
## Read data with 'data.table::fread'
input <- fread("wu03ew_v1.csv", select = 1:3)
setnames(input, 1:3, new = c("origin", "destination","total"))
## Coordinates
centroids <- fread("msoa_popweightedcentroids.csv")
## 'Code' is the key to be used in the joins
setkey(centroids, Code)
## Key of centroids matches `origin` in `input`
origin <- centroids[input[,.(origin, total)]]
setnames(origin, c('East', 'North'), c('xstart', 'ystart'))
## Key of centroids matches `destination` in `input`
destination <- centroids[input[,.(destination)]]
setnames(destination, c('East', 'North'), c('xend', 'yend'))
## Bind both results
trajects <- cbind(origin, destination)
### GRAPHICS SECTION
library(lattice)
library(classInt)
## Background set to black
myTheme <- simpleTheme()
myTheme$background$col <- 'black'
## Palette and classes
nClasses <- 5
pal <- colorRampPalette(c('gray70', 'white'))(nClasses)
classes <- classIntervals(trajects[total > 10, total],
n = nClasses, style = 'quantile')
classes <- findCols(classes)
xyplot(North ~ East, data = centroids,
pch = '.', col = 'lightgray',
aspect = 'iso',
par.settings = myTheme,
panel = function(…){
## panel.xyplot displays the 'centroids'
panel.xyplot(…)
## panel.segments displays the lines using a `data.table`
## query.
trajects[total > 10,
panel.segments(xstart, ystart, xend, yend,
col = pal[classes],
alpha = 0.05, lwd = 0.3)
]
})

view raw

mappingFlows.R

hosted with ❤ by GitHub

This is the result:

london

Related articles
  • Spatial data visualization with R
  • Maps with R (I)
  • Maps with R (II)
  • Maps with R (III)
  • Stamen maps with spplot
Anuncio publicitario

Comparte / Share

  • Correo electrónico
  • Imprimir
  • Twitter
  • Facebook
  • LinkedIn
  • Pinterest
  • Más
  • Tumblr
  • Reddit

Me gusta esto:

Me gusta Cargando...

Label placement with spplot and lattice

29 sábado Dic 2012

Posted by Oscar Perpiñán Lamigueiro in R-english, R-project, visualization

≈ 6 comentarios

Etiquetas

label, lattice, lines, maptools, points, R, spplot

The package maptools includes new functions to label points and labels.

Line labelling

The lineLabel function produces and draws text grobs following the paths defined by a list of Line objects. The sp.lineLabel methods use this function to work easily with spplot. Let’s use the meuse data to illustrate how to use it:

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


data(meuse)
coordinates(meuse) = ~x+y
data(meuse.riv)
## Create a SpatialLines object with ID=1
meuse.sl <- SpatialLines(list(Lines(list(Line(meuse.riv)), "1")))

library(RColorBrewer)
myCols <- adjustcolor(colorRampPalette(brewer.pal(n=9, 'Reds'))(100), .85)

## label is a wrapper to build a character vector with names
## according to the SpatialLines ID's
labs <- label(meuse.sl, 'Meuse River')

With textloc you can choose where the label must be located. For example, with maxDepth the label will start at the point where the line reaches its maximum depth:

sl1 <- list('sp.lineLabel', meuse.sl, label=labs,
            position='below', textloc='maxDepth',
            spar=.2,
            col='darkblue', cex=1,
            fontfamily='Palatino',
            fontface=2)

spplot(meuse.grid["dist"],
       col.regions=myCols, 
       sp.layout = sl1)

https://procomun.files.wordpress.com/2012/12/wpid-linelabelmaxdepth.jpg

The default is to place the label in a stable region where the sign of the slope remains constant.

sl2 <- modifyList(sl1, list(textloc = 'constantSlope'))

spplot(meuse.grid["dist"],
       col.regions=myCols, 
       sp.layout = sl2)

https://procomun.files.wordpress.com/2012/12/wpid-linelabelconstantslope.jpg

But you can define the label location with a numeric index relative to the line length:

sl3 <- modifyList(sl1, list(textloc = 140, position='above'))

spplot(meuse.grid["dist"],
       col.regions=myCols, 
       sp.layout = sl3)

https://procomun.files.wordpress.com/2012/12/wpid-linelabelnumeric.jpg

There is a more sophisticated example in the Spatial Data chapter of my forthcoming book.

Point labelling

This package already provided the pointLabel function with optimization routines to find good locations for point labels without overlaps. This function, useful for base graphics, has been adapted to work with lattice graphics and with the spplot functions.

You will find the panel.pointLabel to use with xyplot and the rest of the lattice family:

library(maptools)
library(lattice)

n <- 15
x <- rnorm(n)*10
y <- rnorm(n)*10
labels <- as.character(round(x, 5))


myTheme <- list(add.text=list(
                  cex=0.7,
                  col='midnightblue',
                  fontface=2,
                  fontfamily='mono'))

xyplot(y~x,
       labels=labels,
       par.settings=myTheme, 
       panel=function(x, y, labels, ...){
         panel.xyplot(x, y, ...)
         panel.pointLabel(x, y, labels=labels, ...)
       })

https://procomun.files.wordpress.com/2012/12/wpid-pointlabel.jpg

And there is the sp.pointLabel method to be combined with spplot:

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

library(RColorBrewer)
myCols <- adjustcolor(colorRampPalette(brewer.pal(n=9, 'Reds'))(100), .85)

pts <- spsample(meuse.grid, n=15, type="random")

## Let's print the name of some R authors 
Rauthors <- readLines(file.path(R.home("doc"), "AUTHORS"))[9:28]
someAuthors <- Rauthors[seq_along(pts)]

sl1 <- list('sp.points', pts, pch=19, cex=.8, col='midnightblue')
sl2 <- list('sp.pointLabel', pts, label=someAuthors,
            cex=0.7, col='midnightblue',
            fontfamily='Palatino')

spplot(meuse.grid["dist"], col.regions=myCols, sp.layout=list(sl1, sl2))

https://procomun.files.wordpress.com/2012/12/wpid-sppointlabel.jpg

There is a more sophisticated example in the Spatial Data chapter of my forthcoming book.

Comparte / Share

  • Correo electrónico
  • Imprimir
  • Twitter
  • Facebook
  • LinkedIn
  • Pinterest
  • Más
  • Tumblr
  • Reddit

Me gusta esto:

Me gusta Cargando...

Spatial Data visualization with R

27 martes Nov 2012

Posted by Oscar Perpiñán Lamigueiro in R-english, R-project, visualization

≈ 1 comentario

Etiquetas

lattice, maps, open data, openstreetmap, R, spatial data, tooltip, visualization

I have published the first version of the code and main figures of the «Spatial Data» chapter of the forthcoming book «Displaying time series, spatial and space-time data with R«.

In this chapter I have included a collection of different techniques. For example, you will find a proportional symbol map with tooltips (click on this image to open the map with tooltips)

or a reference map combining data from OpenStreetMap and hill shading.

Click me!

Comparte / Share

  • Correo electrónico
  • Imprimir
  • Twitter
  • Facebook
  • LinkedIn
  • Pinterest
  • Más
  • Tumblr
  • Reddit

Me gusta esto:

Me gusta Cargando...

Maps with R (II)

20 lunes Feb 2012

Posted by Oscar Perpiñán Lamigueiro in R-english, visualization

≈ 7 comentarios

Etiquetas

GeoTIFF, land cover, lattice, population, R, raster, rasterVis

In my my last post I described how to produce a multivariate choropleth map with R. Now I will show how to create a map from raster files. One of them is a factor which will group the values of the other one. Thus, once again, I will superpose several groups in the same map.

NOTE: Although the procedure described in this post is valid, there is a newer code version in one of the chapters of the book «Displaying time series, spatial and space-time data with R«

First let’s load the packages.

library(raster)
library(rasterVis)
library(colorspace)

Now, I define the geographical extent to be analyzed (approximately India and China).

ext <- extent(65, 135, 5, 55)

The first raster file is the population density in our planet, available at this NEO-NASA webpage (choose the Geo-TIFF floating option, ~25Mb). After reading the data with raster I subset the geographical extent and replace the 99999 with NA.

pop <- raster('875430rgb-167772161.0.FLOAT.TIFF')
pop <- crop(pop, ext)
pop[pop==99999] <- NA
pTotal <- levelplot(pop, zscaleLog=10, par.settings=BTCTheme)
pTotal

The second raster file is the land cover classification (available at this NEO-NASA webpage)

landClass <- raster('241243rgb-167772161.0.TIFF')
landClass <- crop(landClass, ext)

The codes of the classification are described here. In summary, the sea is labeled with 0, forests with 1 to 5, shrublands, grasslands and wetlands with 6 to 11, agriculture and urban lands with 12 to 14, and snow and barren with 15 and 16. This four groups (sea is replaced NA) will be the levels of the factor. (I am not sure if these sets of different land covers is sensible: comments from experts are welcome!)

EDIT: Following a question from a user of rasterVis I include some lines of code to display this qualitative variable in the map.
EDIT2: raster and rasterVis are able to work with categorical data using ratify.

  landClass[landClass %in% c(0, 254)] <- NA
  landClass <- cut(landClass, c(0, 5, 11, 14, 16))
  ## Add a Raster Atribute Table and define the raster as categorical data
  landClass <- ratify(landClass)
  ## Configure the RAT: first create a RAT data.frame using the
  ## levels method; second, set the values for each class (to be
  ## used by levelplot); third, assign this RAT to the raster
  ## using again levels
  rat <- levels(landClass)[[1]]
  rat$classes <- c('Forest', 'Land', 'Urban', 'Snow')
  levels(landClass) <- rat

  levelplot(landClass, col.regions=terrain_hcl(4))

This histogram shows the distribution of the population density in each land class.

s <- stack(pop, landClass)
layerNames(s) <- c('pop', 'landClass')
histogram(~log10(pop)|landClass, data=s,
            scales=list(relation='free'))

Everything is ready for the map. I will create a list of trellis objects with four elements (one for each level of the factor). Each of these objects is the representation of the population density in a particular land class. I use the same scale for all of them to allow for comparisons (the at argument of levelplot receives the correspondent at values from the global map)

at <- pTotal$legend$bottom$args$key$at

pList <- lapply(1:nClasses, function(i){
  landSub <- landClass
  landSub[!(landClass==i)] <- NA
  popSub <- mask(pop, landSub)
  step <- 360/nClasses
  pal <- rev(sequential_hcl(16, h = (30 + step*(i-1))%%360))
  pClass <- levelplot(popSub, zscaleLog=10, at=at,
                      col.regions=pal, margin=FALSE)
})

And that’s all. The rest of the code is exactly the same as in the previous post. If you execute it you will get this image (click on it for higher resolution).

Related articles
  • Maps with R (I)
  • Maps with R (III)
  • Spatial data visualization with R
  • Stamen maps with spplot

https://github.com/oscarperpinan/spacetime-vis/tree/master/code/raster.R

Comparte / Share

  • Correo electrónico
  • Imprimir
  • Twitter
  • Facebook
  • LinkedIn
  • Pinterest
  • Más
  • Tumblr
  • Reddit

Me gusta esto:

Me gusta Cargando...

Maps with R (I)

18 sábado Feb 2012

Posted by Oscar Perpiñán Lamigueiro in R-english, visualization

≈ 58 comentarios

Etiquetas

choropleth, elections, lattice, maps, New York Times, R

This is the first post of a short series to show some code I have learnt to produce maps with R.

NOTE: Although the procedure described in this post is valid, there is a newer code version in one of the chapters of the book «Displaying time series, spatial and space-time data with R«

Some time ago I found this infographic from The New York Times (via this page) and I wondered how a multivariate choropleth map could be produced with R. Here is the code I have arranged to show the results of the last Spanish general elections in a similar fashion.

Some packages are needed:

library(maps)
library(maptools)
## EDITED: if you have rgeos installed you won't need
## gpclibPermit() below.
library(sp)
library(lattice)
library(latticeExtra)
library(colorspace)

Let’s start with the data, which is available here (thanks to Emilio Torres, who «massaged» the original dataset, available from here).

Each region of the map will represent the percentage of votes obtained by the predominant political option. Besides, only four groups will be considered: the two main parties («PP» and «PSOE»), the abstention results («ABS»), and the rest of parties («OTH»).

load(url('http://dl.dropbox.com/u/40293713/spainVotes/votos2011.rda'))
votesData <- votos2011[, 12:1023] ##I don't need all the columns

votesData$ABS <- with(votos2011, Total.censo.electoral - Votos.validos) ##abstention

Max <- apply(votesData, 1, max)
whichMax <- apply(votesData,  1, function(x)names(votesData)[which.max(x)])

## OTH for everything but PP, PSOE and ABS
whichMax[!(whichMax %in% c('PP',  'PSOE', 'ABS'))] <- 'OTH'

## Finally, I calculate the percentage of votes with the electoral census
pcMax <- Max/votos2011$Total.censo.electoral * 100

The Spanish administrative boundaries are available as shapefiles at the INE webpage (~70Mb):

espMap <- readShapePoly(fn="mapas_completo_municipal/esp_muni_0109")
Encoding(levels(espMap$NOMBRE)) <- "latin1"

##There are some repeated polygons which can be dissolved with
## unionSpatialPolygons.
## EDITED: gpclibPermit() is needed for unionSpatialPolygons to work
## but can be ommited if you have rgeos installed
## (recommended, see comment of Roger Bivand below).
espPols <- unionSpatialPolygons(espMap, espMap$PROVMUN)

(EDITED, following the question of Sandra). Spanish maps are commonly displayed with the Canarian islands next to the peninsula. First we have to extract the polygons of the islands and the polygons of the peninsula.

canarias <-  sapply(espPols@polygons, function(x)substr(x@ID, 1, 2) %in% c("35",  "38"))
peninsulaPols <- espPols[!canarias]
islandPols <- espPols[canarias]

Then we shift the coordinates of the islands:

dy <- bbox(peninsulaPols)[2,1] - bbox(islandPols)[2,1]
dx <- bbox(peninsulaPols)[1,2] - bbox(islandPols)[1,2]

islandPols2 <- elide(islandPols, shift=c(dx, dy))
bbIslands <- bbox(islandPols2)

and finally construct a new object binding the shifted islands with the peninsula:

espPols <- rbind(peninsulaPols, islandPols2)

The last step before drawing the map is to link the data with the polygons:

IDs <- sapply(espPols@polygons, function(x)x@ID)
idx <- match(IDs, votos2011$PROVMUN)

##Places without information
idxNA <- which(is.na(idx))

##Information to be added to the SpatialPolygons object
dat2add <- data.frame(prov = votos2011$PROV,
poblacion = votos2011$Nombre.de.Municipio,
Max = Max,  pcMax = pcMax,  who = whichMax)[idx, ]

row.names(dat2add) <- IDs
espMapVotes <- SpatialPolygonsDataFrame(espPols, dat2add)

## Drop those places without information
espMapVotes <- espMapVotes[-idxNA, ]

So let’s draw the map. I will produce a list of plots, one for each group. The «+.trellis» method of the latticeExtra package with Reduce superposes the elements of this list and produce a trellis object. I will use a set of sequential palettes from the colorspace package with a different hue for each group.

classes <- levels(factor(whichMax))
nClasses <- length(classes)

pList <- lapply(1:nClasses, function(i){
  mapClass <- espMapVotes[espMapVotes$who==classes[i],]
  step <- 360/nClasses ## distance between hues
  pal <- rev(sequential_hcl(16, h = (30 + step*(i-1))%%360)) ## hues equally spaced
  pClass <- spplot(mapClass['pcMax'], col.regions=pal, lwd=0.1,
                   at = c(0, 20, 40, 60, 80, 100))
  })

p <- Reduce('+', pList)

However, the legend of this trellis object is not valid.

First, a title for the legend of each element pList will be
useful. Unfortunately, the levelplot function (the engine under the
spplot method) does not allow for a title with its colorkey
argument. The frameGrob and packGrob of the grid package will do the work.

addTitle <- function(legend, title){
  titleGrob <- textGrob(title, gp=gpar(fontsize=8), hjust=1, vjust=1)
  legendGrob <- eval(as.call(c(as.symbol(legend$fun), legend$args)))
  ly <- grid.layout(ncol=1, nrow=2, widths=unit(0.9, 'grobwidth', data=legendGrob))
  fg <- frameGrob(ly, name=paste('legendTitle', title, sep='_'))
  pg <- packGrob(fg, titleGrob, row=2)
  pg <- packGrob(pg, legendGrob, row=1)
  }

for (i in seq_along(classes)){
  lg <- pList[[i]]$legend$right
  lg$args$key$labels$cex=ifelse(i==nClasses, 0.8, 0) ##only the last legend needs labels
  pList[[i]]$legend$right <- list(fun='addTitle',
                                  args=list(legend=lg, title=classes[i]))
}

Now, every component of pList includes a legend with a title below. The last step is to modify the legend of the p trellis object in order to merge the legends from every component of pList.

## list of legends
legendList <- lapply(pList, function(x){
  lg <- x$legend$right
  clKey <- eval(as.call(c(as.symbol(lg$fun), lg$args)))
  clKey
})

##function to pack the list of legends in a unique legend
##adapted from latticeExtra::: mergedTrellisLegendGrob
packLegend <- function(legendList){
  N <- length(legendList)
  ly <- grid.layout(nrow = 1,  ncol = N)
  g <- frameGrob(layout = ly, name = "mergedLegend")
  for (i in 1:N) g <- packGrob(g, legendList[[i]], col = i)
  g
}

## The legend of p will include all the legends
p$legend$right <- list(fun = 'packLegend',  args = list(legendList = legendList))

Here is the result with the provinces boundaries superposed (only for the peninsula due to a problem with the definition of boundaries the Canarian islands in the file) and a rectangle to separate the Canarian islands from the
rest of the map (click on the image to get a SVG file):

provinces <- readShapePoly(fn="mapas_completo_municipal/spain_provinces_ag_2")

canarias <- provinces$PROV %in% c(35, 38)
peninsulaLines <- provinces[!canarias,]

p +
  layer(sp.polygons(peninsulaLines,  lwd = 0.5)) +
  layer(grid.rect(x=bbIslands[1,1], y=bbIslands[2,1],
                  width=diff(bbIslands[1,]),
                  height=diff(bbIslands[2,]),
                  default.units='native', just=c('left', 'bottom'),
                  gp=gpar(lwd=0.5, fill='transparent')))

https://procomun.files.wordpress.com/2012/04/spainvotes.jpeg

Related articles
  • Maps with R (II)
  • Maps with R (III)
  • Spatial data visualization with R
  • Stamen maps with spplot

Comparte / Share

  • Correo electrónico
  • Imprimir
  • Twitter
  • Facebook
  • LinkedIn
  • Pinterest
  • Más
  • Tumblr
  • Reddit

Me gusta esto:

Me gusta Cargando...
← Entradas anteriores

Omnia sunt communia!

“As we enjoy great advantages from the inventions of others, we should be glad of an opportunity to serve others by any invention of ours, and this we should do freely and generously”

RSS Feed RSS - Entradas

Introduce tu correo electrónico para suscribirte a este blog y recibir avisos de nuevas entradas.

Únete a otros 91 suscriptores

Posts Más Vistos

  • Libro de Energía Solar Fotovoltaica
  • Vector fields with streamlines
View Oscar Perpiñán Lamigueiro's profile on LinkedIn
profile for Oscar Perpiñán on Stack Exchange, a network of free, community-driven Q&A sites
Citations for Oscar Perpinan

@oscarperpinan

Mis tuits

Enlaces

  • Displaying time series, spatial and space-time data with R
  • Introducción a R

Top Clicks

  • github.com/oscarperpinan/…

Etiquetas

actuar anarquía autocontención ayuda mutua bloch ciencia circuitos eléctricos común cooperación CRAN creative commons data data.table data analysis ecología emacs Energía energía solar fotovoltaica entropía EOI esperanza exergía forecast fotovoltaica galeano github GNU General Public License información investigación jorge riechmann kropotkin latex lattice lewis mumford libro lyx map maps memoir meteorological data movimientos sociales nestoria nicolas georgescu roegen NWP open data openstreetmap orgmode pensar photovoltaics pobreza R R-project radiación solar raster rasterVis reproducible research siar small multiples software libre solar solar radiation sp spatial spatial data spplot tecnología time series trellis tufte uned utopía Vector field visualización visualization WRF

RSS R-bloggers

  • Simulating confounders, colliders and mediators by @ellis2013nz
  • Introduction to Mixed-effects Models in R workshop
  • Learn to ‘Make an outstanding Shiny App’ with us
  • Sorting, Ordering, and Ranking: Unraveling R’s Powerful Functions
  • Galton and Watson voluntarily skipping some generations
  • {attachment} v0.4.0: Breaking changes and configuration file for a better experience
  • Sharing the Big Book of R upgrade proposal
  • The do.call() function in R: Unlocking Efficiency and Flexibility
  • Version 1.0.0 of NIMBLE released, providing automatic differentiation, Laplace approximation, and HMC sampling
  • Checking normality in R

RSS dataanalytics

  • Se ha producido un error; es probable que la fuente esté fuera de servicio. Vuelve a intentarlo más tarde.

RSS Solar Energy (Elsevier)

  • Se ha producido un error; es probable que la fuente esté fuera de servicio. Vuelve a intentarlo más tarde.

RSS Progress in Photovoltaics

  • Se ha producido un error; es probable que la fuente esté fuera de servicio. Vuelve a intentarlo más tarde.

RSS Madridmasd

  • Cambio Climático, Gobiernos y Empresas; Casi todas sus promesas son mentira
  • Proyectos estratégicos para la recuperación y transformación económica (PERTE)
  • ¿En quién confiar? Debate sobre IA, ética y comunicación en el congreso de medios de comunicación
  • La Civilización perdida bajo los suelos y vegetación de Latinoamérica.
  • IMDEA Software e IMDEA Networks trabajan para desplegar en la Comunidad de Madrid “MadQCI”: la mayor red cuántica de Europa
  • Etnografía y Paleoecología de los Paisajes Precolombinos: El Papel del Fuego en las Culturas Aborígenes del SO de América del Norte (tribus Karuk y Yurok)
  • Homenaje a Manuel Espadas Burgos
  • Super Mario World
  • Taller doctoral en Historia de las Relaciones Internacionales 2023 en la Universidad de Padua
  • Directiva Europea de Protección de Suelos: Dos décadas de “vistas” y “consideraciones” y otras desconsideraciones
Licencia Creative Commons
Salvo indicación en contra todos los contenidos están publicados por Oscar Perpiñán Lamigueiro bajo una Licencia Creative Commons Reconocimiento-No Comercial-Compartir Igual 3.0 Unported.

Blog de WordPress.com.

Privacidad y cookies: este sitio utiliza cookies. Al continuar utilizando esta web, aceptas su uso.
Para obtener más información, incluido cómo controlar las cookies, consulta aquí: Política de cookies
  • Seguir Siguiendo
    • Omnia sunt Communia!
    • Únete a 91 seguidores más
    • ¿Ya tienes una cuenta de WordPress.com? Inicia sesión.
    • Omnia sunt Communia!
    • Personalizar
    • Seguir Siguiendo
    • Regístrate
    • Iniciar sesión
    • Denunciar este contenido
    • Ver sitio web en el Lector
    • Gestionar las suscripciones
    • Contraer esta barra
 

Cargando comentarios...
 

    A %d blogueros les gusta esto: