R: Como obter latitudes e longitudes de um RasterLayer?

14

Eu sou um iniciante absoluto de dados geográficos; portanto, perdoe-me se a pergunta não for apropriada.

Baixei dados do NCDC NARR e consegui carregar no R usando o rasterpacote. Gostaria de obter uma lista com latitude, longitude e valor. Entendo que rasterToPoints()deve fazer exatamente o que quero, no entanto, meus valores de latitude e longitude parecem estranhos:

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

Suponho que devo fazer algo com a projeção que atualmente é Lambert Conformal Conic (LCC). Aqui estão mais informações sobre a varredura.

> 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

O que devo fazer para obter valores reais de latitude e longitude nos EUA?

janosdivenyi
fonte

Respostas:

14

você realmente precisa reprojetar a varredura em uma projeção geográfica (graus decimais) usando "projectRaster" ou "spTransform". Observe também as definições de CRS sp que especificam a sequência de projeção desejada. O exemplo da ajuda para o "projectRaster" é bastante claro em como fazer isso.

Se você coagir seus dados raster em um objeto SpatialPointsDataFrame, use "spTransform" e puxe as coordenadas do slot @coordinates e as adicione ao data.frame no slot @data. Aqui está um exemplo de como isso seria.

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)

Devo observar que não é uma boa prática converter rasters em uma classe de objeto de vetor e nega as vantagens do pacote raster, fornecendo um processamento seguro de memória. Muitas vezes, é prudente pensar realmente no seu problema e avaliar se você o está abordando corretamente. Se o OP tivesse fornecido o contexto do motivo pelo qual eles precisam de coordenadas [x, y] para cada célula, a comunidade do fórum pode ter sido capaz de fornecer alternativas computacionais que manteriam o problema em um ambiente de varredura.

Jeffrey Evans
fonte
1
Uma maneira de tomar sua precaução (sobre como evitar a conversão de dados) é desprojetar a varredura original (talvez para uma grade muito grossa), criar duas grades de latitude e longitude cobrindo a extensão da desprojeção e projetá-las de volta ao extensão da grade original. Nenhuma classe de vetor é criada: é inteiramente um conjunto de operações de varredura.
whuber
4

Obtenha as coordenadas dos centros celulares e crie um objeto espacial:

spts <- rasterToPoints(r, spatial = TRUE)

Transforme os pontos no seu alvo desejado:

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

Os valores já são copiados como colunas neste SpatialPointsDataFrame.

print(llpts)

Agora, para finalizar, obtenha um data.frame:

x <- as.data.frame(llpts)

Há uma implementação geral disso no pacote SGAT, veja a função lonlatFromCellaqui:

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

mdsumner
fonte
Eu tentei isso, mas recebi a seguinte mensagem de erro: > 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
Na verdade, você não precisa transferir os atributos, eu removo isso.
Mdsumner #
Além do conselho do pacote SGAT, essa não é exatamente a mesma resposta / exemplo que eu forneci? As coordenadas não são propagadas para o data.frame no slot de dados, apenas os valores da varredura. As coordenadas são, de fato, mantidas no slot de coordenadas e precisam ser adicionadas ao data.frame.
10135 Jeffrey Evans
Obrigado, adicionei a etapa as.data.frame. Eu acho que é um péssimo conselho adicionar as coordenadas como atributos - especialmente movendo o slot - já que as coordenadas do objeto podem mudar. Se você quiser um data.frame bruto, basta criar um. Eu não me importo onde está a informação, talvez apenas edite a sua e podemos zap esta resposta.
Mdsumner #
O OP queria especificamente coordenadas e acho que é redundante salvar em um data.frame separado. Normalmente, não gosto de adicionar coordenadas ao slot de dados principalmente porque é redundante com o slot de coordenadas. Fora isso, não é "péssimo conselho" adicionar informações ao slot de dados. E se você quiser ter dois sistemas de coordenadas. Você pode adicionar lat / long ao slot de dados e ter o objeto em uma projeção totalmente diferente. Além disso, se você deseja apenas exportar um arquivo simples, e não um formato GIS, por si só, é possível adicionar coordenadas ao data.frame e salvar como um CSV.
Jeffrey Evans
0

Parece que você tem coordenadas projetadas lá (não as coordenadas de latitude / longitude, também conhecidas como GCS). Provavelmente não estava claro para você que esse era o problema. Veja este post. Convertendo sistema de coordenadas geográficas em R

jbchurchill
fonte
Não peguei o post que você mencionou antes de responder. Você pode sinalizar isso como duplicado. Embora a adição da coerção SpatialPointsDataFrame e da atribuição de coordenadas seja um pouco diferente. Sua chamada.
Jeffrey Evans
Pensei em marcá-lo como tal, mas pensei que se outra pessoa estiver procurando por uma resposta semelhante sem saber que precisa projetar os valores, isso pode aparecer para eles. Além de sua resposta oferecer uma maneira diferente de chegar lá (votado).
precisa saber é o seguinte
Eu tentei olhar para as fontes que você listou. A fim de obter latitudes / longitudes padrão que eu emitei lonlat_r <- projectRaster(r, crs="+init=epsg:4326"). No entanto, a extensão da nova varredura está -181.3232, 181.4938, -1.590457, 87.76154 (xmin, xmax, ymin, ymax)longe do que eu esperaria dos EUA (que deveria estar entre 30 a 70 e -60 a -160). Eu deveria ter entendido algo errado.
precisa saber é o seguinte