Проекция объектов пузыря / круга на картах

Question

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

https://www.ngdc.noaa.gov/nndc/struts/results?type_0=Exact&query_0=$ID&t=101650&s=13&d=189&dfn=signif.txt

Затем запустил следующий код для очистки / организации этого набора данных:

library(dplyr)
library(tmap)
library(sf)
earthquake<-read.table("signif.txt",sep="\t",header=TRUE,fill=TRUE) %>% filter(!is.na(LATITUDE) & !is.na(LONGITUDE)) %>% st_as_sf(coords=c("LONGITUDE","LATITUDE"))

Затем запустил следующеекод для отображения карты всех землетрясений магнитудой 9 и более:

tmap_mode("view")
tm_shape(earthquake %>% filter(EQ_PRIMARY > 9))+tm_bubbles(size = "EQ_PRIMARY",col="red",popup.vars=c("EQ_PRIMARY"))

Я получаю это сообщение об ошибке, поскольку я никогда не назначал проекцию на данные: Currect projection of shape earthquake %>% filter(EQ_PRIMARY > 9) unknown. Long-lat (WGS84) is assumed. Это нормально, и я получаю прикрепленную картинку:

Проблема заключается в том, что величина этого землетрясения на Аляске на самом деле составляет 9,2, тогда как сила южного в Чили составляет 9,5,все же круг Аляски заметно больше!Значки пузырьков дальше от экватора проецируются и искажаются под проекцией Меркатора.

Поэтому я пытаюсь изменить проекцию моих данных на LAEA:

st_crs(earthquake)<-"+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs "

Но теперь, когда я бегута же карта, что и выше, круги отображаются в правильном размере, но базовая карта не отображается, так как я полагаю, что tmap не имеет базовой карты LAEA?Вот где я заблудился.

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

Какое здесь решение?Я хотел бы использовать красивую карту Меркатора, потому что она выглядит красиво, но я не хочу, чтобы такие вещи, как символы, искажались ею.Нужно ли определять новый столбец размера, чтобы противодействовать искажению Меркатора, например earthquake %>% mutate(EQ_PRIMARY1 = EQ_PRIMARY / (abs(LATITUDE)+1)), но заменить его фактической исследуемой функцией, которая будет противодействовать эффекту размера?Это распространенная проблема в этой области, или этот пакет просто не работает правильно?

Answer 1

Я нашел некоторый обходной путь, как я описал, получив масштабный коэффициент для каждой точки на основе ее широты.Простите, если это странный обходной путь, но я немного новичок в этом.Базовая процедура была такой:

1. Добавьте две новые геометрии, основанные на проекции существующей геометрии - спроецируйте одну из них на EQC, а одну - ту же самую проекцию EQC, но сместив вверх на 1 единицу (я сделал это путемустановка ложного севера на 1 в proj4).В основном нам нужны две геометрии, где одна единица находится к северу от другой.
2. Преобразуйте эти две геометрии в веб-меркатор, используя st_transform.
3. Измерьте расстояние по оси y между получающимися геометриями.Это то, насколько Web Mercator масштабирует 1 единицу расстояния на этой конкретной широте.Это очень близко к 1.000 для точек, близких к экватору, и почти к 2.5 около полюсов.
4. Масштабирование размера точек с помощью обратного квадрата вышеуказанного фактора.

Кодвыглядит следующим образом:

earthquake$yeqc0<-earthquake %>% st_transform("+proj=eqc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs") %>% st_geometry()
earthquake$yeqc1<-earthquake %>% st_transform("+proj=eqc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=1 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs") %>% st_geometry()

Обратите внимание на тонкое y_0 = 1 выше.

st_crs(earthquake$yeqc1)<-st_crs(earthquake$yeqc0)
earthquake$ymerc<-st_distance(st_transform(earthquake$yeqc0,crs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"),st_transform(earthquake$yeqc1,crs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"),by_element=TRUE)
earthquake<-earthquake %>% mutate(EQ_PRIMARY1 = EQ_PRIMARY * (1 / ymerc)^2)
tm_shape(earthquake %>% filter(EQ_PRIMARY > 9))+tm_bubbles(size = "EQ_PRIMARY1",col="red",popup.vars=c("EQ_PRIMARY","EQ_PRIMARY1"))

Эта карта создаетимеет точки как правильный размер.Код немного неуклюжий, но может быть и хуже.Возможно, что-то, что заслуживает пользовательской функции, если ее не существует.

Answer 2

Во-первых, вам нужно установить crs для проекции, в которой он находится в данный момент, вероятно, WGS84. Затем используйте st_transform, чтобы изменить координаты. Вот код, но он, похоже, показывает аналогичный результат. Возможно, сработает другая проекция, или попробуйте изменить проекцию в функции tm_shape.

library(dplyr)
library(tmap)
library(sf)
earthquake<-read.table("https://www.ngdc.noaa.gov/nndc/struts/results?type_0=Exact&query_0=$ID&t=101650&s=13&d=189&dfn=signif.txt",
                       sep="\t",header=TRUE,fill=TRUE) %>%
  filter(!is.na(LATITUDE) & !is.na(LONGITUDE)) %>%
  st_as_sf(coords=c("LONGITUDE","LATITUDE"))
st_crs(earthquake) <- 4326
earthquake <- st_transform(earthquake, "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")
tmap_mode("view")
tm_shape(earthquake %>% filter(EQ_PRIMARY > 9))+
  tm_bubbles(size = "EQ_PRIMARY",col="red",popup.vars=c("EQ_PRIMARY"))