7

I'm trying to plot a map of the Pacific using World2Hires in R's mapproj library but there's an odd glitch when I try to fill the countries. How do I fix this?

library(maps)
library(mapproj)
library(mapdata)
map("world2Hires", 
    xlim=c(120, 260), 
    ylim=c(-60, 40), 
    boundary=TRUE, 
    interior=TRUE,
    fill=TRUE,
    col="gray30",
)
map.axes()

Here's the output:

Broken Map Image

Puzzled79
  • 1,037
  • 2
  • 10
  • 23
  • For others looking for a solution to this issue, there are couple of other answers of interest. For a solution that leaves the polygons intact, allowing a global map centered on the Pac, see http://stackoverflow.com/a/10749877/3897439 from Josh O'Brien. To split the polygons near the dateline, see http://stackoverflow.com/a/5538551/3897439 by Joris Meys. – Cotton.Rockwood Aug 04 '14 at 16:52

2 Answers2

7

The issue seems to be with a small subset of areas which cause wrapping. From some trial and error saving the original map call like mapnames <- map(...) and then passing subsets of this list to the regions= argument in a new call, I could avoid the wrapping around of the fills. E.g.:

library(maps)
library(mapproj)
library(mapdata)
map("world2Hires", regions=mapnames$names[c(1:7,14:641)],
    xlim=c(120, 260), 
    ylim=c(-60, 40), 
    boundary=TRUE, 
    interior=TRUE,
    fill=TRUE
)
map.axes()

enter image description here

As to a more thorough or sensible solution to prevent this happening, I am stumped. Playing with the wrap= option does nothing helpful, and likewise for the other options. As a side-note, this issue does not appear using the "world" database, but only pops up for "world2" and "world2Hires".

thelatemail
  • 91,185
  • 12
  • 128
  • 188
  • @Puzzled79 - apparently this issue has been discussed previously here: http://stackoverflow.com/questions/5353184/fixing-maps-library-data-for-pacific-centred-0-360-longitude-display – thelatemail Feb 10 '14 at 05:36
5

The answer by @thelatemail is the best and simplest solution I've seen to this problem. To make it more universal, though, it is better to remove the polygons by name. This is because depending on the limits you give your first call to map(), the indices of the polygon names can be different.

library(maps)
library(mapproj)
library(mapdata)

mapnames <- map("world2Hires", xlim=c(120, 260), ylim=c(-60, 40),
                fill=TRUE, plot=FALSE)

mapnames2 <- map("world2Hires", xlim=c(100, 200), ylim=c(-20, 60),
                 fill=TRUE, plot=FALSE)

mapnames$names[10]
[1] "Mali"
mapnames2$names[10]
[1] "Thailand"

There are 8 countries that are intersected by the Prime Meridian: United Kingdom, France, Spain, Algeria, Mali, Burkina Faso, Ghana, and Togo. By matching these country names with mapnames$names, you can remove the polygons regardless of your original extent:

remove <- c("UK:Great Britain", "France", "Spain", "Algeria", "Mali",
            "Burkina Faso", "Ghana", "Togo")

map("world2Hires", regions=mapnames$names[!(mapnames$names %in% remove)],
    xlim=c(120, 260), 
    ylim=c(-60, 40), 
    boundary=TRUE, 
    interior=TRUE,
    fill=TRUE
)
map.axes()

You could also use grepl() but because polygons are named heirarchically, you may remove some sub-polygons of the nations in question. For example, mapnames$names[grepl("UK", mapnames$names)] returns 34 matches.

I would have suggested this as an edit, but I don't have privileges yet.

Cotton.Rockwood
  • 1,601
  • 12
  • 29
  • This is a great addition expanding on the previous answer. Identifying the specific subset of problematic countries is quite helpful. If you keep posting useful things like this, your lack of privileges won't be for long. – thelatemail Aug 01 '14 at 00:31
  • Thanks, @thelatemail! It also occurred to me that there could be a solution wherever your map break is. This might be done by overlaying a N-S line (at your break longitude) on your polygons to identify which polygons it intersects. Perhaps over(lines, polygons)? – Cotton.Rockwood Aug 01 '14 at 04:07
  • There is this answer http://stackoverflow.com/questions/5353184/fixing-maps-library-data-for-pacific-centred-0%C2%B0-360%C2%B0-longitude-display which rearranges polygons that split across the longitude line. – thelatemail Aug 01 '14 at 04:08
  • I think @Joris Meys answer that you linked above is still probably the most elegant solution, but not as simple as yours. Another interesting solution is from kohske: http://stackoverflow.com/questions/10620862/use-different-center-than-the-prime-meridian-in-plotting-a-world-map. – Cotton.Rockwood Aug 01 '14 at 04:14
  • It does not work for the whole world (xlim=c(0, 360), ylim=c(-90, 90)) – user3910073 Apr 28 '16 at 10:23