1. Pergunta

Encontrei um aviso usando a função projectRaster () no pacote raster em R. Um exemplo completo e reproduzível é colado abaixo.

   Warning message:
   In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1],  :
   33940 projected point(s) not finite

Minha pergunta é: Esse aviso é um problema que preciso corrigir se estiver trabalhando com dados terrestres? Em outras palavras, os dados são "perdidos". Este seria um grande problema para mim, se for. Se for esse o caso, você sabe se existe uma maneira de corrigi-lo?

Eu procurei uma solução para este problema online e encontrei alguma menção a ele aqui , aqui e aqui, mas acho que nenhum oferece uma resposta adequada para esse problema.

2. Carregue a biblioteca raster

  library(raster)

3. Primeiro, crie um mapa (longo lat) do mundo, com uma constante em cada célula da grade

Estou colocando uma constante em cada célula da grade para ver se consigo diagnosticar o problema de onde o aviso provavelmente afetará os dados, se houver.

 rastertest.longlat<-raster(ncol=360, nrow=180)
 values(rastertest.longlat)<-c(rep(1,n=180*360))

4. Se você reprojetar o mapa (longo lat) em uma grade de área igual, ele produzirá a mensagem de aviso

   rastertest.eck4<-projectRaster(rastertest.longlat, res=c(100000,100000), crs="+proj=eck4",method="ngb", over=T)

Warning message:
 In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1],  :
33940 projected point(s) not finite

Penso que esta mensagem está basicamente dizendo que a re-projeção falhou para algumas das células da grade.

5. Mas se você traçar os dois mapas, não parecerá que esse aviso cause problemas para os dados

Ou seja, você não vê lacunas em branco nos dados plotados. Meu palpite é que as células perdidas são células não terrestres, no topo e na borda do mundo. Alguma ideia?

par(mfrow=c(2,1))

plot(rastertest.longlat, col="blue")
data(wrld_simpl)
wrld <- spTransform(wrld_simpl, CRS('+proj=longlat'))
plot(wrld, add=TRUE)

plot(rastertest.eck4, col="blue") 
wrld <- spTransform(wrld_simpl, CRS('+proj=eck4'))
plot(wrld, add=TRUE)

Resposta curta: está tudo bem, e você pode agradecer rasterpor continuar a concluir, em vez de falhar, e por informar que alguns dados foram perdidos.

Resposta longa:

isso dependerá da projeção e, neste caso, provavelmente está apenas nas "bordas". qual é a aresta e como ela se manifesta para uma determinada instância de uma determinada família de projeção é a parte "depende".

Você pode ver que não são os pontos centrais das células perdidas:

tpoints <- rgdal::project(coordinates(rastertest.longlat), "+proj=eck4")
sum(is.na(tpoints))
#[1] 0

Mas provavelmente são os cantos e, possivelmente, as bordas da própria célula. Talvez isso mostre que projetos rasterizados se baseiam na extensão das células, não apenas em seus pontos centrais.

 rgdal::project(as.matrix(expand.grid(x = c(-180, 0, 180), y = c(-90,0, 90))), "+proj=eck4")

Admito que esperava que fosse de onde vêm os valores ausentes, então talvez projectRasteresteja se estendendo um pouco mais ao norte e ao sul? Defina valores lá para latitude fora do intervalo -90/90 e você começará a receber o aviso. Vou acompanhar se tiver a chance de explorar mais.

Finalmente, você provavelmente deve usar um parâmetro elipsóide ou de dado explícito, ou seja, "+ proj = eck4 + ellps = WGS84".

Esta não é uma resposta completa à minha pergunta original, em termos de todos os detalhes, mas fornece ao leitor interessado um lugar para ir.

Em geral, os dados são perdidos e distorcidos durante a reprojeção de rasters de projeções longlat para áreas iguais. Isso pode ser problemático para a análise global. Se você pode obter seus dados em formato vetorial, é melhor reprojetar usando esse formato.

Aqui está uma referência sobre o problema geral. E aqui outra sobre tentar quantificar a perda. Espero que ajude.