Agregando pontos à grade usando R

14

Eu tenho uma pergunta em relação à agregação espacial em R. O que estou tentando fazer é agregar um conjunto de dados de pontos a uma grade. No entanto, não tenho certeza de como fazer isso, pois tenho pouca experiência com esse tipo de coisa. Eu esperava que algum de vocês tivesse alguma orientação útil / uma possível solução.

Meu ponto de vista é um conjunto de dados que contém dados georreferenciados sobre eventos de conflito na África (consulte www.acleddata.com). Os pontos são georreferenciados com coordenadas de latitude / longitude e contêm dados sobre o tipo e a hora do evento. O que eu quero fazer é agregar esses pontos em uma grade de 1x1 grau.

Assim, uma célula de grade deve conter as informações dos pontos de dados se um evento ocorrer dentro dessa célula de grade. O produto final disso deve ser um quadro de dados ou algo que eu possa exportar para um arquivo csv, pois os dados devem ser usados ​​em um conjunto de dados de painel para análise estatística.

Até agora, carreguei e plotei os dados e o shapefile usando o código abaixo. Acredito que devo usar a função over do pacote sp para agregar, mas não sei como. Espero que um de vocês possa ajudar.

O código que eu usei até agora pode ser encontrado aqui com o resultado visual correspondente ali .

Sugestões para fazer isso no QGIS também são bem-vindas.

horseoftheyear
fonte
Esta é uma operação simples e rápida que requer nada mais que um pouco de aritmética. Mas em qual formato você deseja a saída? "CSV" sugere apenas que deve ser uma tabela relacional, mas isso apresenta um problema: quando você agrega, cada célula corresponde potencialmente a um número variável de pontos. Geralmente, você seleciona uma das duas opções: você gera um registro por ponto (incluindo o ID da célula que o contém) ou um registro por célula e inclui alguns resumos estatísticos dos pontos que ele contém. De que você precisa?
whuber
1
Desculpe, eu não especifiquei isso. O que eu preciso é de um registro por célula . Eu uso o arquivo csv para criar dados do painel no formato de ano da célula .
precisa saber é o seguinte

Respostas:

12

Os dados baixados contêm alguns erros de localização francos; portanto, a primeira coisa a fazer é limitar as coordenadas a valores razoáveis:

data.df <- read.csv("f:/temp/All_Africa_1997-2011.csv", header=TRUE, sep=",",row.names=NULL)
data.df <- subset(data.df, subset=(LONGITUDE >= -180 & LATITUDE >= -90))

O cálculo das coordenadas e identificadores das células da grade é apenas uma questão de truncar os decimais dos valores de latitude e longitude. (Em geral, para rasters arbitrários, primeiro centralize e dimensione-os para unidade de tamanho de célula, trunque as casas decimais e depois redimensione e atualize novamente a posição original, conforme mostrado no código jiabaixo.) Podemos combinar essas coordenadas em identificadores exclusivos, anexando-os ao quadro de dados de entrada e grave o quadro de dados aumentado como um arquivo CSV. Haverá um registro por ponto:

ji <- function(xy, origin=c(0,0), cellsize=c(1,1)) {
  t(apply(xy, 1, function(z) cellsize/2+origin+cellsize*(floor((z - origin)/cellsize))))
}
JI <- ji(cbind(data.df$LONGITUDE, data.df$LATITUDE))
data.df$X <- JI[, 1]
data.df$Y <- JI[, 2]
data.df$Cell <- paste(data.df$X, data.df$Y)

Você pode desejar uma saída que resuma eventos dentro de cada célula da grade. Para ilustrar isso, vamos calcular as contagens por célula e produzir essas, um registro por célula:

counts <- by(data.df, data.df$Cell, function(d) c(d$X[1], d$Y[1], nrow(d)))
counts.m <- matrix(unlist(counts), nrow=3)
rownames(counts.m) <- c("X", "Y", "Count")
write.csv(as.data.frame(t(counts.m)), "f:/temp/grid.csv")

Para outros resumos, altere o functionargumento no cálculo de counts. (Como alternativa, use o software de planilha ou banco de dados para resumir o primeiro arquivo de saída pelo identificador de célula.)

Como verificação, vamos mapear as contagens usando os centros de grade para localizar os símbolos do mapa. (Os pontos localizados no Mar Mediterrâneo, na Europa e no Oceano Atlântico têm locais suspeitos: suspeito que muitos deles resultem da mistura de latitude e longitude no processo de entrada de dados.)

count.max <- max(counts.m["Count",])
colors = sapply(counts.m["Count",], function(n) hsv(sqrt(n/count.max), .7, .7, .5))
plot(counts.m["X",] + 1/2, counts.m["Y",] + 1/2, cex=sqrt(counts.m["Count",]/100),
     pch = 19, col=colors,
     xlab="Longitude of cell center", ylab="Latitude of cell center",
     main="Event counts within one-degree grid cells")

Mapa da áfrica

Este fluxo de trabalho está agora

  • Completamente documentado (por meio do Rpróprio código),

  • Reproduzível (executando novamente este código),

  • Extensível (modificando o código de maneiras óbvias) e

  • Razoavelmente rápido (toda a operação leva menos de 10 segundos para processar essas observações 53052).

whuber
fonte
O código é perfeitamente reproduzível. Eu tenho uma pergunta adicional embora. Em vez de um resumo, como anexar as informações do arquivo de dados de entrada à célula na grade criada?
precisa saber é o seguinte
1
Isso não é possível com uma tabela de saída , porque as informações completas para células têm comprimento variável. A maneira correta de registrar isso é com a primeira forma de saída que exibi: um registro por ponto com um atributo identificador de célula. Um desses dois formatos - as tabelas por ponto e por célula - será esperado por qualquer programa estatístico que você esteja usando.
whuber
1
Ah ok. Eu vejo o que você quer dizer. Só é necessário criar uma grade para todas as células e mesclá-la. Obrigado pela ajuda.
precisa saber é o seguinte
3

Bem, o que você quer é um chamado "Junção Espacial" básico, que combine dois arquivos de forma e aloque a soma (número da contagem) à tabela de atributos resultante. Se você procurar "Junção espacial no R", encontrará vários exemplos, mesmo aqui no GIS.Stackexchange. Eu rapidamente pesquisei e encontrei, por exemplo, este código publicado em uma lista de discussão.

Se você deseja obter uma junção de atributo espacial no QGIS, faça o seguinte:

  • Salve suas formas como arquivos .shp (comando writeOGR do pacote rgdal)
  • Carregue-os no QGIS. Recrie sua grade de vetores por meio do plug-in MMQGIS (Criar -> Criar camada de grade) com a escala apropriada.
  • Use a ferramenta "Join Attributes" no menu Vector -> Data Management. Selecione um atributo da sua camada de pontos (pode ser uma coluna simples que representa valores TRUE (1) ou FALSE (0) para diferentes eventos de conflito).
  • Selecione sua grade e Soma todas as ocorrências e execute. Depois, eu também cortaria sua grade com uma forma do continente africano.

Se o Join de alguma forma falhar (não funciona para mim sempre), siga o SEXTANTE e procure a caixa de ferramentas SAGA, que também possui funções de junção muito boas.

Maçarico
fonte
Embora essa seja uma solução, é particularmente complexa e ineficiente, uma vez que resumir os pontos de uma grade é apenas uma questão de algumas operações aritméticas simples, que Rse destacam em. O uso de shapefiles, rgdalQGIS e Sextante é como recomendar que alguém contrate uma moderna planta industrial automatizada para unir duas placas :-).
whuber
Vou tentar essa abordagem neste fim de semana. Num futuro próximo, talvez eu queira combinar vários arquivos de forma entre si, para que isso possa ser útil. Obrigado pela contribuição e sugestões.
precisa saber é o seguinte
@ whuber: Isso é verdade, mas se você deseja distribuir e talvez estilizar sua saída, um shapefile é a escolha óbvia. No entanto, bom exemplo de R!
Curlew
Eu finalmente tentei. Mas o problema dessa abordagem é que ela soma todas as observações ao polígono. Embora, idealmente, queira manter as informações sobre diferentes eventos ao longo do tempo. Mas pode ser que eu tenha feito algo errado.
horseoftheyear