1

I've got a map in ggplot (made up of geom_polygon, lines, points and segments). The original was in WGS, with all units in degrees lat/lon. All of the additional data (for lines and points) are also in degrees. I want to retransform the entire thing into an equal-area projection (Mollweide). The first step (transforming the underlying rworldmap polygon) was fairly painless, although there are some odd horizontal lines running through. However, I'm not sure how to tackle the next step. As I said, all the data are in degrees, and I want to specify xlim/ylim, axes labels and (curved) gridlines in degrees. Will I have to convert everything into Mollweide meters, or is it somehow possible to just reshape my final map object in one painless step?

This is what i have so far, to demonstrate these weird horizontal lines.

library(rworldmap)
library(rgdal)
library(ggplot2)

moll_crs<-CRS("+proj=moll +R=10567000 +lon_0=0 +x_0=0 +y_0=0 +units=m +towgs84=0,0,0,0,0,0,0 +no_defs")

ggplot() + 
  geom_polygon(data = spTransform(getMap(resolution = 'low'), CRSobj = moll_crs), 
               aes(x = long, 
                   y = lat, 
                   group = group),
               fill = 'gray90', 
               colour = 'gray10', 
               size = 0.3) +
  coord_fixed()

This results in horizontal lines.

enter image description here

EDIT:

Following @hrbmstr 's answer, I'm using his GeoJSON file with coord_map("molleweide"):

ggplot() + 
  geom_map(data=world, map=world,
           aes(x=long, y=lat, map_id=id), 
           fill="gray90", color="gray10", size=0.3) +
  coord_map("mollweide",
            xlim = c(-50, 40))

The xlim argument is messing things up though, adding in horizontal lines. I thought the coord_***(xlim) arguments just changed the viewable area, but it seems to be affecting how the map is drawn. Any ideas? enter image description here

Bhargav Rao
  • 50,140
  • 28
  • 121
  • 140
TFinch
  • 361
  • 2
  • 9

2 Answers2

3

You can let ggplot do it for you, and you can also get rid of the lines by using a better map. The GeoJSON file I reference below was created from the Natural Earth shapefiles and you can just use the shapefile from there (or the optimized GeoJSON I made from here).

To show how it deals with additional points, I added close centroids for USA & AUS

library(sp)
library(ggplot2)
library(rgdal)
library(rgeos)

world <- readOGR("ne_50m_admin_0_countries.geojson", "OGRGeoJSON")
outline <- bbox(world)
outline <- data.frame(xmin=outline["x","min"],
                      xmax=outline["x","max"],
                      ymin=outline["y","min"],
                      ymax=outline["y","max"])

world <- fortify(world)

points <- data.frame(lon=c(-98.35, 134.21), lat=c(39.5, -25.36))

gg <- ggplot()
gg <- gg + geom_rect(data=outline, 
                     aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), 
                     color=1, fill="white", size=0.3)
gg <- gg + geom_map(data=world, map=world,
                    aes(x=long, y=lat, map_id=id), 
                    fill="gray90", color="gray10", size=0.3)
gg <- gg + geom_point(data=points, aes(x=lon, y=lat), size=3)
gg <- gg + labs(x=NULL, y=NULL)
gg <- gg + coord_map("mollweide")
gg <- gg + theme_bw()
gg <- gg + theme(panel.grid=element_blank())
gg <- gg + theme(panel.border=element_blank())
gg <- gg + theme(axis.ticks=element_blank())
gg <- gg + theme(axis.text=element_blank())
gg


gg

enter image description here

You can clip the world ahead of time to avoid the issues with coord_map:

world <- readOGR("ne_50m_admin_0_countries.geojson", "OGRGeoJSON")
clipper <- as(extent(-50, 40, -60, 60), "SpatialPolygons")
proj4string(clipper) <- CRS(proj4string(world))

world <- gIntersection(world, clipper, byid=TRUE)
world <- fortify(world)

gg <- ggplot()
gg <- gg + geom_map(data=world, map=world,
                    aes(x=long, y=lat, map_id=id), 
                    fill="gray90", color="gray10", size=0.3)
gg <- gg + labs(x=NULL, y=NULL)
gg <- gg + coord_map("mollweide")
gg

(NOTE: I didn't include all the code from above to avoid clutter). You can also use the same bounding box code to put a projected border around the clipped map, too.

enter image description here

hrbrmstr
  • 77,368
  • 11
  • 139
  • 205
  • Note you should never use a projection like this without a background showing where the edge of the world is. Cartographers will cry. – Spacedman Jan 12 '15 at 12:28
  • OP is going to use `xlim`/`ylim` to zoom in somewhere, so I didn't think it was important for the example, but you're right. will fix. – hrbrmstr Jan 12 '15 at 12:35
  • Thanks for this - very neat. There's an issue with using `xlim` (inside `coord_map`) though. If the limits add up to 0 (e.g. `xlim = c(-40, 40)` it works fine, but with `xlim = c(-50, 40)` I'm getting horizontal lines again (using @hrbrmstr 's GeoJSON). No problems with `ylim` – TFinch Jan 12 '15 at 12:47
  • If you crop first - http://stackoverflow.com/questions/13982773/crop-for-spatialpolygonsdataframe/13986029#13986029 - you should be fine – hrbrmstr Jan 12 '15 at 12:52
1

Excluding antarctica from the rworldmap map is one way of correcting the annoying horizontal lines. Plus note that you only need the coarse resolution map which you get from the default value of resolution in rworldmap::getMap().

Here is a solution to the first part of your problem with minimal changes to your code.

moll_crs<-CRS("+proj=moll +ellps=WGS84")
#first get countries excluding Antarctica which can crash spTransform
#for a world plot you only need coarse resolution which is the default for getMap()
sPDF <- getMap()[getMap()$ADMIN!='Antarctica',]

ggplot() + 
  geom_polygon(data = spTransform(sPDF, CRSobj = moll_crs), 
               aes(x = long, 
                   y = lat, 
                   group = group),
               fill = 'gray90', 
               colour = 'gray10', 
               size = 0.3) +
  coord_fixed()
Andy
  • 1,821
  • 13
  • 23