使用墨卡托投影在ggplot2中设置地图限制

与此
question非常相关我试图绘制一些世界区域,现在使用墨卡托投影,但在添加x和y限制时一直遇到麻烦:

ggplot(world, mapping = aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "black", colour = "black") +
  coord_map(projection = "mercator", xlim = c(-125, -30), ylim = c(-60, 35))

《使用墨卡托投影在ggplot2中设置地图限制》
显然不是很好.当我使用coord_cartesian(如建议的here)设置限制时,我松开了墨卡托投影:

ggplot(world, mapping = aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "black", colour = "black") +
  coord_map(projection = "mercator") +
  coord_cartesian(xlim = c(-125, -30), ylim = c(-60, 35))

《使用墨卡托投影在ggplot2中设置地图限制》

当我使用lims时,我得到了我想要的拉丁美洲:

ggplot(world, mapping = aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "black", colour = "black") +
  coord_map(projection = "mercator") +
  lims(x = c(-125, -30), y = c(-60, 35))

《使用墨卡托投影在ggplot2中设置地图限制》

问题是,这种方法并不总是有效,例如在非洲或亚洲,我开始对情节边界采取一些疯狂行为:

ggplot(world, mapping = aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "black", colour = "black") +
  coord_map(projection = "mercator") +
  lims(x = c(-20, 45), y = c(-50, 40))
  # lims(x = c(40, 150), y = c(-10, 55))

《使用墨卡托投影在ggplot2中设置地图限制》
《使用墨卡托投影在ggplot2中设置地图限制》

最佳答案 这里的解决方案可以是将lat / lon坐标转换为“正确的”web mercator坐标(这里我使用的是epsg 3857,这是“google”投影),然后使用这些“新”坐标进行绘图.

假设原始坐标是latlon wgs84(epsg 4326),可以这样实现:

worldmerc <-  SpatialPointsDataFrame(coords = data_frame(x = world$long, y = world$lat), 
                                data = world, proj4string = CRS("+proj=longlat +datum=WGS84")) %>%
         subset((lat < 90 & lat > -90)) %>%   # needed because transform not defined at the poles !!!!
         spTransform(CRS("+init=epsg:3857")) 
worldmerc  <-  mutate(worldmerc@data, longmerc = coordinates(worldmerc)[,1], latmerc = coordinates(worldmerc)[,2])

绘制整个数据给出了这一点(注意使用coord_fixed来保持纵横比!:

ggplot(worldmerc, mapping = aes(x = longmerc, y = latmerc, group = group)) +
  geom_polygon(fill = "black", colour = "black") +coord_fixed()

《使用墨卡托投影在ggplot2中设置地图限制》

现在,问题是要进行子集化,你现在需要输入“map”坐标而不是lat long,但也可以调整:

#For South America
xlim = c(-125, -30)
ylim = c(-60, 35)

lims = SpatialPoints(coords = data_frame(x = xlim, y = ylim), proj4string = CRS("+proj=longlat +datum=WGS84"))%>%
  spTransform(CRS("+init=epsg:3857"))

ggplot(worldmerc, mapping = aes(x = longmerc, y = latmerc, group = group)) +
  geom_polygon(fill = "black", colour = "black")+
  coord_fixed(xlim = coordinates(lims)[,1], ylim = coordinates(lims)[,2])

《使用墨卡托投影在ggplot2中设置地图限制》

#for africa
xlim = c(-20,45)
ylim = c(-50,40)

lims = SpatialPoints(coords = data_frame(x = xlim, y = ylim), proj4string = CRS("+proj=longlat +datum=WGS84"))%>%
       spTransform(CRS("+init=epsg:3857"))

ggplot(worldmerc, mapping = aes(x = longmerc, y = latmerc, group = group)) +
  geom_polygon(fill = "black", colour = "black")+
  coord_fixed(xlim = coordinates(lims)[,1], ylim = coordinates(lims)[,2])

《使用墨卡托投影在ggplot2中设置地图限制》

如您所见,在这两种情况下,您都可以获得“正确”的地图.

现在,您可能要做的最后一件事可能是在轴上有“lat / lon”坐标.这有点像黑客但可以这样做:

library(magrittr)
xlim = c(-125, -30)
ylim = c(-60, 35)

# Get the coordinates of the limits in mercator projection
lims = SpatialPoints(coords = data_frame(x = xlim, y = ylim), 
                     proj4string = CRS("+proj=longlat +datum=WGS84"))%>%
       spTransform(CRS("+init=epsg:3857"))

# Create regular "grids" of latlon coordinates and find points 
# within xlim/ylim - will be our labels

majgrid_wid_lat = 20
majgrid_wid_lon = 30

majbreaks_lon = data_frame(x=seq(-180,  180, majgrid_wid_lon)) %>% 
                filter(x >= xlim[1] & x <= xlim[2]) %>% 
                as.data.frame()
majbreaks_lat = data_frame(x=seq(-90,   90, majgrid_wid_lat)) %>%
                filter(x >= ylim[1] & x <= ylim[2]) %>% 
                as.data.frame()

#Find corresponding mercator coordinates

mercbreaks_lat = SpatialPoints(coords = expand.grid(x = majbreaks_lon$x, y = majbreaks_lat$x), proj4string = CRS("+init=epsg:4326"))%>%
  spTransform(CRS("+init=epsg:3857")) %>% coordinates() %>% extract(,2) %>% unique() 
mercbreaks_lon = SpatialPoints(coords = expand.grid(x = majbreaks_lon$x, y = majbreaks_lat$x), proj4string = CRS("+init=epsg:4326"))%>%
  spTransform(CRS("+init=epsg:3857")) %>% coordinates()  %>% extract(,1) %>% unique()

# Plot using mercator coordinates, but latlon labels

ggplot(worldmerc, mapping = aes(x = longmerc, y = latmerc, group = group)) +
  geom_polygon(fill = "black", colour = "black") +
  coord_fixed(xlim = coordinates(lims)[,1], ylim = coordinates(lims)[,2])+
  scale_x_continuous("lon", breaks = mercbreaks_lon, labels = signif(majbreaks_lon$x, 2)) + 
  scale_y_continuous("lat", breaks = mercbreaks_lat, labels = signif(majbreaks_lat$x,2))+theme_bw() 

, 这使:

《使用墨卡托投影在ggplot2中设置地图限制》

它有点复杂,可能有更好的方法,但它可以解决问题,并且可以很容易地在函数中进行转换.

HTH,

洛伦佐

点赞