Построение широтных / длинных координат в карте растрового слоя (фактора) формального класса с проекцией aea в R?

#r #spatial #r-raster

#r #пространственный #r-растр

Вопрос:

Я довольно новичок в работе с пространственными фреймами данных, и у меня есть то, что я считал относительно простой задачей: возьмите фрейм данных из 6 точек, со столбцами x и y, представляющими позиции широты / длины этих точек, и спроецируйте их так, чтобы их можно было использовать в пространственном фрейме данных, который у меня естьсделано.

Вот как я закодировал в 6 точках:

 d1 <- structure(list(latitude = c(37.427733, 37.565759, 37.580956, 37.429285, 37.424270, 37.502496), longitude = c(-108.011061, -107.814039, -107.676662, -107.677166, -108.898826, -108.586042)))

d2 <- as.data.frame(d1)

d3 <- SpatialPointsDataFrame(c(d2[,c('longitude','latitude')]), data = d2)

 

И я попытался изменить / назначить проекцию для них (эти данные широты / длины были взяты из Google Maps), но, похоже, я не могу заставить это работать. Проекция для данных, на которые я хочу наложить эти точки, выглядит следующим образом:

  proj=aea  lat_0=23  lon_0=-96  lat_1=29.5  lat_2=45.5  x_0=0  y_0=0  datum=WGS84  units=m  no_defs
 

Итак, в основном мой вопрос: как я могу преобразовать эти широты / длины в формат для x / y, который использует эта проекция? Вот объем набора данных, который я хочу наложить на него для справки, показывая, что он явно не в простой широте / длине:

 class      : Extent 
xmin       : -1145835 
xmax       : -1011345 
ymin       : 1613205 
ymax       : 1704855 
 

Заранее большое вам всем спасибо!

Ответ №1:

Вам нужно (повторно) спроецировать свои пространственные точки в ту же проекцию, что и ваш другой источник данных. Я больше знаком с sf пакетом для работы с пространственной векторной информацией в R, но похоже, что вы используете sp packcage.

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

 library(sp)

d1 <- structure(list(latitude = c(37.427733, 37.565759, 37.580956, 37.429285, 37.424270, 37.502496), longitude = c(-108.011061, -107.814039, -107.676662, -107.677166, -108.898826, -108.586042)))

d2 <- as.data.frame(d1)

d3 <- SpatialPointsDataFrame(c(d2[,c('longitude','latitude')]), data = d2)


proj4string(d3) <- CRS(" init=epsg:4326")

sf::st_bbox(d3)
# xmin       ymin       xmax       ymax 
# -108.89883   37.42427 -107.67666   37.58096 
 

Глядя на сводку или этот экстент, вы можете видеть, что широта и длина соответствуют оригиналу. Теперь мы можем перепроектировать, используя предоставленную вами proj4string

 target_crs = CRS(" proj=aea  lat_0=23  lon_0=-96  lat_1=29.5  lat_2=45.5  x_0=0  y_0=0  datum=WGS84  units=m  no_defs")  # This is your string assigned to an object


d4 = spTransform(d3, target_crs) ## object used to transform your data frame

sf::st_bbox(d4)
# xmin     ymin     xmax     ymax 
# -1127246  1661667 -1018854  1679942 
 

Ваши координаты теперь представлены в целевом пространстве и должны быть нанесены поверх ваших фоновых данных.

Если вы хотите использовать альтернативный пакет, sf::st_transform() это немного проще в использовании, и sf пакет, как правило, более удобен для пользователя.

Комментарии:

1. Большое вам спасибо! Это отлично работает. Я думаю, вы правы, что формат sp несколько сбил меня с толку. Большое спасибо за добрые и полезные отзывы.