R: Как получить широты и долготы от RasterLayer?

14

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

Я загрузил данные из NCDC NARR и сумел загрузить их в R с помощью rasterпакета. Я хотел бы получить список с широтой, долготой и значением. Я понимаю, что rasterToPoints()должен делать именно то, что я хочу, однако мои значения широты и долготы выглядят странно:

r <- raster(myfile)
data_matrix <- rasterToPoints(r)
head(data_matrix)
            x       y value
[1,] -5405401 4347242    70
[2,] -5372938 4347242    88
[3,] -5340475 4347242    76
[4,] -5308012 4347242    85
[5,] -5275549 4347242    87
[6,] -5243086 4347242    88

Я полагаю, что я должен сделать что-то с проекцией, которая в настоящее время является конформной конической формой Ламберта (LCC). Вот дополнительная информация о растре.

> r
class       : RasterLayer 
dimensions  : 277, 349, 96673  (nrow, ncol, ncell)
resolution  : 32463, 32463  (x, y)
extent      : -5648874, 5680713, -4628777, 4363474  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=50 +lat_2=50 +lat_0=50 +lon_0=-107 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs 
data source : mypath-to-file
names       : value

Что мне делать, чтобы получить реальные значения широты и долготы США?

janosdivenyi
источник

Ответы:

14

вам нужно на самом деле перепроецировать растр в географическую проекцию (в десятичных градусах), используя «projectRaster» или «spTransform». Также посмотрите на определения CRS sp, которые определяют желаемую строку проекции. Пример в справке для «projectRaster» совершенно ясно, как это сделать.

Если вы приведете свои растровые данные в объект SpatialPointsDataFrame, то вы будете использовать «spTransform» и извлекать координаты из слота @coordinates и добавлять их в data.frame в слоте @data. Вот пример того, как это будет выглядеть.

library(raster)
library(rgdal) # for spTransform

# Create data
r <- raster(ncols=100, nrows=100)
  r[] <- runif(ncell(r))
  crs(r) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
  projection(r)

# Convert raster to SpatialPointsDataFrame
r.pts <- rasterToPoints(r, spatial=TRUE)
  proj4string(r.pts)

# reproject sp object
geo.prj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" 
r.pts <- spTransform(r.pts, CRS(geo.prj)) 
  proj4string(r.pts)

# Assign coordinates to @data slot, display first 6 rows of data.frame
r.pts@data <- data.frame(r.pts@data, long=coordinates(r.pts)[,1],
                         lat=coordinates(r.pts)[,2])                         
head(r.pts@data)

Следует отметить, что преобразование растров в класс векторных объектов нецелесообразно и сводит на нет преимущества пакета растров, обеспечивающего безопасную обработку памяти. Часто разумно подумать о своей проблеме и оценить, правильно ли вы к ней подходите. Если ОП предоставил контекст, объясняющий, почему им нужны координаты [x, y] для каждой ячейки, форумное сообщество могло бы предоставить вычислительные альтернативы, которые позволили бы сохранить проблему в растровой среде.

Джеффри Эванс
источник
1
Один из способов проявить осторожность (избегая преобразования данных) состоит в том, чтобы отменить проецирование исходного растра (возможно, на очень грубую сетку), создать две сетки значений широты и долготы, охватывающие степень непроекции, и спроецировать их обратно в степень исходной сетки. Векторные классы не создаются: это полностью набор растровых операций.
whuber
4

Получить координаты центров ячеек и создать пространственный объект:

spts <- rasterToPoints(r, spatial = TRUE)

Преобразуйте точки к желаемой цели:

library(rgdal)
llprj <-  "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
llpts <- spTransform(spts, CRS(llprj))

Значения уже скопированы как столбцы в этом SpatialPointsDataFrame.

print(llpts)

Теперь, чтобы закончить, получите data.frame:

x <- as.data.frame(llpts)

Есть общая реализация этого в пакете SGAT, смотрите lonlatFromCellздесь функцию :

https://github.com/SWotherspoon/SGAT/blob/master/R/Raster.R

mdsumner
источник
Я попробовал это, но получил следующее сообщение об ошибке: > llpts$layer1 <- values(r[[1]]) Error in [[<-. Data.frame (* tmp *, name, value = c(NA, NA, NA, NA, NA, : replacement has 96673 rows, data has 95025
janosdivenyi
На самом деле вам не нужно передавать атрибуты, я удалю это.
mdsumner
Кроме рекомендаций пакета SGAT, это не тот самый ответ / пример, который я предоставил? Координаты не передаются в data.frame в слоте данных, только значения из растра. По сути, координаты хранятся в слоте координат и должны быть добавлены в data.frame.
Джеффри Эванс
Спасибо, я добавил шаг as.data.frame. Я думаю, что это ужасный совет, чтобы добавить координаты в качестве атрибутов - особенно с помощью слота - поскольку координаты объекта могут измениться. Если вам нужен необработанный data.frame, просто создайте его. Мне все равно, где информация, может быть, просто отредактируйте вашу, и мы можем убрать этот ответ.
mdsumner
ОП специально хотел координаты, и я думаю, что это избыточно для сохранения в отдельном data.frame. Обычно я не люблю добавлять координаты в слот данных в основном, потому что он избыточен со слотом координат. Кроме этого, это не «ужасный совет» для добавления информации в слот данных. Что делать, если вы хотите иметь две системы координат. Вы можете добавить lat / long в слот данных и получить объект в совершенно другой проекции. Кроме того, если вы хотите просто экспортировать плоский файл, а не формат ГИС как таковой, вы можете добавить координаты в data.frame и сохранить как CSV.
Джеффри Эванс
0

Похоже, что у вас там есть проекционные координаты (а не координаты широты / долготы или координаты GCS) Возможно, вам не было ясно, в чем проблема. Смотрите этот пост. Преобразование географической системы координат в R

jbchurchill
источник
Я не уловил пост, на который вы ссылаетесь, прежде чем я ответил. Вы можете пометить это как дубликат. Хотя добавление приведения и координат в SpatialPointsDataFrame немного отличается. Ваш звонок.
Джеффри Эванс
Я думал о том, чтобы пометить его как таковой, но подумал, что если другой человек ищет аналогичный ответ, не зная, что ему нужно спроецировать значения, это может показаться ему. Кроме того, ваш ответ предлагает другой способ добраться (Upvoted).
jbchurchill
Я пытался посмотреть на источники, которые вы перечислили. Для того чтобы получить стандартные широты / долготы я выдал lonlat_r <- projectRaster(r, crs="+init=epsg:4326"). Тем не менее, размер нового растра -181.3232, 181.4938, -1.590457, 87.76154 (xmin, xmax, ymin, ymax)далек от того, что я ожидал бы от США (который должен быть где-то между 30 до 70 и от -60 до -160). Я должен был что-то неправильно понять.
Janosdivenyi