Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Plot a geographical map with symbols according to frequencies

Tags:

r

gis

qgis

I would like to plot a geographical map with n numbers of squares according to the frequency of a disease within a county. Like this one here: enter image description here

But I couldn't figure out how to do this with R or qGIS. Thanks for y'all help.

like image 364
Dominik Avatar asked Dec 16 '22 06:12

Dominik


2 Answers

First write a little function called 'stackbox' that plots (using "rect") the squares for the stack in the right place.

Here's the first line of that function:

stackbox <- function(x,y,n,size,maxheight=5){
  • where 'size' is the height and width of the boxes, and 'maxheight' lets you have a stack that is 5 high or 10 high or whatever.

Then call that function for every county that has cases.

Where exactly were you stuck in the process?

Here's the function in full:

stackbox <- function(x,y,n,size,maxheight=5,...){
  stackheight = seq(0,n,by=maxheight)
  stackheight=diff(unique(c(stackheight,n)))

  for(col in 1:length(stackheight)){
    xl=rep(x+(col-1)*size,stackheight[col]) - (length(stackheight)/2)*size
    yb=y+size*((1:stackheight[col])-1) - (max(stackheight)/2)*size
    xr=xl+size
    yt=yb+size
    rect(xl,yb,xr,yt,...)
  }
}

Example:

plot(1:10)
for(i in 1:10){
 stackbox(i,i,i,3,size=.1,col="red",border="white")
}

To do this on a map, you need the sp and maptools packages, and a shapefile or other geospatial data source that has your data in it:

africa=readShapeSpatial(file.path(mapLib,"africa.shp"))
plot(africa,border="gray")
coords=coordinates(africa)
for(i in 1:nrow(africa)){
  if(cases[i]>0){
    stackbox(coords[i,1],coords[i,2],africa$cases[i],1,border="#606060",col="#0083FE")
  }
}

Map using stacked boxes

I've picked colours that look a bit like your original. Note that the boxes are in plot coordinates, so I had to make them 0.1 of a degree. You may wish to transform your map to a euclidean projection (using spTransform from package:gdal) and then they'll be in those units and properly square.

Doing nice text labels like in your original would be tricky though...

like image 126
Spacedman Avatar answered Jan 31 '23 08:01

Spacedman


As an alternative I would suggest a bubble plot, which is quite easy to make using ggplot2 using the point geometry. An example form the help pages (http://had.co.nz/ggplot2/geom_point.html):

p <- ggplot(mtcars, aes(wt, mpg)) 
p + geom_point(aes(size = qsec)) 

which leads to:

ggplot bubble plot

The size of the bubble represents the amount measured. The example does not use geographical coordinates, but they can be used easily. The ggplot2 package also supports adding spatial polygons as an additional layer, e.g. using geom_path. The help file of coord_map shows a few examples using map data. See also the fortify function to transform a SpatialPolygons object to a data.frame (which is needed for ggplot2).

like image 28
Paul Hiemstra Avatar answered Jan 31 '23 07:01

Paul Hiemstra