Визуальная ошибка при изменении центрального меридиана проекции Робинсона с помощью ggplot2 - PullRequest
Новая
       63

Визуальная ошибка при изменении центрального меридиана проекции Робинсона с помощью ggplot2

0 голосов
/ 15 мая 2019

Я пытаюсь спроецировать карту мира в проекции Робинсона, где центральный меридиан отличается от 0. Согласно этому потоку StackOverFlow , это должно быть легко (хотя в примере используется sp).

Вот мой воспроизводимый код: 

library(sf)
library(ggplot2)
library(rnaturalearth)

world <- ne_countries(scale = 'small', returnclass = 'sf')

# Notice +lon_0=180 instead of 0
world_robinson <- st_transform(world, crs = '+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

ggplot() +
geom_sf(data = world_robinson)

Это результат.Полигоны замыкаются с одной стороны на другую проекции.

image

Попытка с sp дает тот же эффект.Я также пробовал с шейп-файлом, включающим только полигоны с береговых линий (без политических границ) из http://www.naturalearthdata.com/, и эффект был аналогичным.

image

Я пытался запустить свой фрагментна двух независимых установках R в Mac OS X и Ubuntu 18.04.

Ответы [ 2 ]

6 голосов
/ 15 мая 2019

Это расширение к ответу Z.lin (т.е. сначала используйте этот ответ для вычисления world_robinson).Однако есть еще один полезный шаг, который можно добавить.После проецирования области, состоящие из более чем одного многоугольника, потому что они пересекаются с одной стороны карты на другую в исходной проекции (см. Антарктида, Фиджи и Россия), все еще имеют это разделение после репроекции.Например, вот крупный план Антарктиды, где мы можем видеть, что у нее есть граница на главном меридиане, где ни один из них не должен быть:

enter image description here

КПрошивая эти области вместе, мы можем сначала выяснить, какие полигоны являются проблемами, найдя те, которые пересекают простой меридиан: 

bbox = st_bbox(world_robinson)
bbox[c(1,3)] = c(-1e-5,1e-5)
polygon2 <- st_as_sfc(bbox)

crosses = world_robinson %>%
  st_intersects(polygon2) %>%
  sapply(length) %>%
  as.logical %>%
  which

Теперь мы можем выбрать эти полигоны и установить их размер буфера равным нулю: 

library(magrittr)
world_robinson[crosses,] %<>%
  st_buffer(0) 

ggplot(world_robinson) + geom_sf()

Как мы видим, на карте больше нет расщеплений по главному меридиану:

enter image description here

4 голосов
/ 15 мая 2019

Полигоны, которые пересекают линию меридиана, растягиваются по всей карте после трансформации. Один из способов обойти это - разделить эти многоугольники посередине, чтобы все многоугольники были либо полностью к западу, либо к востоку от линии. 

# define a long & slim polygon that overlaps the meridian line & set its CRS to match
# that of world
polygon <- st_polygon(x = list(rbind(c(-0.0001, 90),
                                     c(0, 90),
                                     c(0, -90),
                                     c(-0.0001, -90),
                                     c(-0.0001, 90)))) %>%
  st_sfc() %>%
  st_set_crs(4326)

# modify world dataset to remove overlapping portions with world's polygons
world2 <- world %>% st_difference(polygon)

# perform transformation on modified version of world dataset
world_robinson <- st_transform(world2, 
                               crs = '+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

# plot
ggplot() +
  geom_sf(data = world_robinson)

plot

...