7

I'm trying to plot points from study sites with a background map of Africa. I can create the two independently, but I am having a hard time overlaying them on top of eachother.

The map of Africa I am using is an Esri shapefile from maplibrary.org. It is available from my dropbox at https://www.dropbox.com/s/etqdw3nky52czv4/Africa%20map.zip. I have the points in a text file, also available from my drop box. https://www.dropbox.com/s/scvymytjsr5pvaf/SPM-437-22Nov12.txt. They refer to studies on molecular drug resistance of malaria parasites. I would like to plot them so that the color is the proportion of parasites with the drug resistant genetic marker and the size is the number of parasites tested.

Plotting the points independently:

qplot(Longitude, Latitude, data = d.spm.437, colour = Frc437, size = Tot437)

Plotting the map of Africa:

library(maptools)
africa = readShapePoly("Africa.shp")
africa.map = fortify(africa, region="COUNTRY")
qplot(long, lat, data = africa.map, geom="path", group=group)

Any help on putting these two together while preserving the display of the points would be appreciated.

David
  • 265
  • 4
  • 13
  • 1
    can you add to your post the eaxact link to the map? – agstudy Dec 01 '12 at 01:02
  • 4
    If you adjust your code to use the `ggplot()` syntax your life will be easier. You can specify separate datasets in your call to `ggplot()`...so your code may look something like `ggplot() + geom_path(data = XX, aes(...)) + geom_point(data = YY, aes(...))` – Chase Dec 01 '12 at 01:08
  • I've been trying to do it that way but am running into problems with the map. For example I tried this code: "ggplot() + geom_path(data = africa.map, eas(long, lat)) + geom_point(d.spm.437, aes(Longitude, Latitude, colour = d.spm.437$Frc437, size = d.spm.437$Tot437))" But iI get the following error: ggplot2 doesn't know how to deal with data of class uneval – David Dec 01 '12 at 01:19
  • `eas(long, lat))`? Not `aes(long, lat))`? – SlowLearner Dec 01 '12 at 01:26
  • That was just a typo. Yes, aes(long, lat) – David Dec 01 '12 at 01:29

1 Answers1

9

Try something like this. Seems to work for me. I think some of your lat-long coordinates are wrong though. The fill colour for geom_point is currently set to Tot437 so you might want to change that.

map

library(ggplot2)
library(rgdal)

africa <- readOGR("c:/test", layer = "Africa")
africa.map = fortify(africa, region="COUNTRY")

africa.points = read.table("c:/test/SPM-437-22Nov12.txt", header = TRUE, sep = ",")
names(africa.points)[which(names(africa.points) == 'Longitude')] <- 'long' # rename lat and long for consistency with shp file
names(africa.points)[which(names(africa.points) == 'Latitude')] <- 'lat'

ggplot(africa.map, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "black", size = 1, fill = "white", aes(group = group)) +
    geom_point(data = africa.points, aes(x = long, y = lat, fill = Tot437, group = NULL), size = 4, shape = 21, colour = "black", size = 3)

Incidentally, looking at your map you may have difficulty getting a good detailed view of individual areas, so one way to tackle that would be by subsetting, in this case with the data frames. You could do this:

africa.map <- africa.map[africa.map$id == 'Madagascar', ]
africa.points <- africa.points[africa.points$Country == 'Madagascar', ]
ggplot(africa.map, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "black", size = 1, fill = "white", aes(group = group)) +
    geom_point(data = africa.points, aes(x = long, y = lat, fill = Tot437, group = NULL), size = 2, shape = 21, colour = "black", size = 2)

...which should get you something similar to this:

madagascar

SlowLearner
  • 7,907
  • 11
  • 49
  • 80
  • nice job; this would be better, though, if you could provide a *relocatable* example (i.e. without any absolute pathnames) ... hence more easily reproducible ... – Ben Bolker Dec 01 '12 at 02:23
  • A fair point, especially for non-Windows people. However, doesn't that imply dumping files into your main `R` folder, thus making rather a mess? – SlowLearner Dec 01 '12 at 02:25
  • Up. I wish there were a way of boosting points for good answers. I know... I'll search your history and bump other good answers (or questions). – IRTFM Dec 01 '12 at 03:31
  • 1
    no, @SlowLearner, it doesn't -- set up a temporary directory wherever you want, put the files there, change the working directory to that location, and go from there. The point is that you change the working directory *outside* of the context of the code you're presenting. (If you want you can show the `setwd()` command as the first line of your code -- that shows everything but still makes it easy for people to skip it if they have their code located somewhere else.) – Ben Bolker Dec 01 '12 at 04:34
  • Thank you SlowLearner for the response. I am getting an error when I try to fortify the shape file though: "Error: TopologyException: found non-noded intersection between LINESTRING (34.7279 1.59723, 34.7278 1.59729) and LINESTRING (34.7278 1.59723, 34.7278 1.59729) at 34.727793021883102 1.5972887049072426" Any idea what that might be? – David Dec 01 '12 at 13:38
  • I figured out a way around the error I was getting - I dissolved the shapefile using QGIS and saved a new shapefile with countries only. Thanks again SlowLearner. – David Dec 01 '12 at 18:58
  • @David Hmmm... That's odd. I used the file you provided with no such problems. OK, well, if the QGIS route works for you then good! – SlowLearner Dec 01 '12 at 19:01