[R] Help with map

[Duncan Murdoch]

On 04/05/2007 9:00 PM, Alberto Vieira Ferreira Monteiro wrote:
> [for those that worry about these things, this _is_ a homework
> assignment. However, it's not an R homework, it's a Geography
> and History homework... and I want to use R to create a pretty
> map]
> 
> Roger Bivand wrote:
>>> Is there any way to associate one color to each country?
>> Try:
>>
>> map_poly_obj <- map("worldHires", c("Argentina", "Brazil"), plot=FALSE,
>>   fill=TRUE)
>> str(map_poly_obj)
>>
>> and you'll see that the component of interest is the named polygons, of
>> which there are 28, namely
>>
> Ok, I guess I can see what you mean. It worked, but I don't think
> this is a practical way to draw things.
> 
> For example, suppose [this would help homework mentioned above] I
> want to draw a series of maps showing the evolution of Communism
> in the XX century. I would like to start with a 1917 map, showing most
> countries as in...
> 
> map("worldHires")
> 
> ... but with the Soviet Union in red. I don't see how I could mix the two 
> maps (BTW, there's no Russia in worldHires, but there is a USSR...)
> 
> map("worldHires"); map("worldHires", "USSR", col="red", fill=T)
> 
>> map_poly_obj$names
>>
>> So you can build a matching colours vector, or:
>>
>> library(sp)
>> library(maptools)
>> IDs <- sapply(strsplit(map_poly_obj$names, ":"), function(x) x[1])
>> SP_AB <- map2SpatialPolygons(map_poly_obj, IDs=IDs,
>>   proj4string=CRS("+proj=longlat +datum=wgs84"))
>>
>> but
>>
>> plot(SP_AB, col=c("cyan", "green"))
>>
>> still misses, because some polygons have their rings of coordinates in
>> counter-clockwise order, so:
>>
>> pl_new <- lapply(slot(SP_AB, "polygons"), checkPolygonsHoles)
>> slot(SP_AB, "polygons") <- pl_new
>> # please forget the assignment to the slot and do not do it unless you can
>> # replace what was there before
>>
>> plot(SP_AB, col=c("cyan", "green"), axes=TRUE)
>>
>> now works. Moreover, SP_AB is a SpatialPolygons object, which can be
>> promoted to a SpatialPolygonsDataFrame object, for a data slot holding a
>> data.frame with row names matching the Polygons ID values:
>>
>> sapply(slot(SP_AB, "polygons"), function(x) slot(x, "ID"))
>>
>> So adding a suitable data frame gets you to the lattice graphics method
>>
>> spplot(SP_AB, "my_var")
>>
>> Hope this helps,
>>
> So, in the above mentioned case, I could do something like:
> 
> library(mapdata)
> commies <- c("USSR", "Mongolia") 
> # Mongolia was the 2nd communist country, in 1925
> map_poly_obj <- map("worldHires", plot=FALSE)
> map_poly_commies <- map("worldHires", commies,
>   plot=FALSE, fill=TRUE)
> plot(map_poly_obj, type="l")
> polygon(map_poly_commies, col="red", border="black")
> 
> I guess I can keep going, unless there is a simpler solution.

Here's the sloppy code I used to put together the map in the persp3d 
example in rgl 0.71.  I've just made a few edits; I hope it still runs.
The idea is to do a lot of calculations on a vector of colours, then
plot them all at once.

> library(mapdata)
> 
> names <- map("worldHires", plot=FALSE, namesonly=TRUE)

> countries <- gsub(":.*", "", names)  # the first part of the name
> lakes <- grep("Lake", names)
> lakes <- lakes[-15]
> lakes <- unique(c(lakes, grep("Ozero", names), 
>                          grep("Vodokh", names),
> 			   grep("Nuur", names),
> 			   grep("Ostrov Ol'khon", names),
> 			   grep("Hayk", names),
> 			   grep("Lago", names)))  # The map doesn't distinguish these
> 
> seas <- grep("Sea", names) 
> nfld <- grep("^Newfoundland$", names) # nfld should be coloured as Canada
> canada <- grep("^Canada$", names)

> svalbard <- grep("Svalbard", names)   # Svalbard coloured as Norway
> norway <- grep("^Norway$", names)
> 
> nums <- as.numeric(factor(countries))
> N <- max(nums)
> set.seed(3)
> col <-  hcl(h=sample((1:N)*360/N),
>             c=sample(seq(35,45,len=N)),
>             l=sample(seq(75,85,len=N)))[nums]  # random colours by country
> col[c(lakes,seas)] <- "white" # with lakes and seas white
> col[nfld] <- col[canada]
> col[svalbard] <- col[norway]
> 
>  png(width=2200, height=2200, file='world.png')
>  map('worldHires', fill = TRUE, col = col, ylim=c(-90,90))
>  abline(h=c(-90, 90))
>  abline(v=c(-180.05, 180.05))
>  dev.off()

I hope this gives you some ideas.

Duncan Murdoch