Reprogrammer le raster en R: avertit que le ou les points projetés ne sont pas finis?

8

1. Question

J'ai rencontré un avertissement en utilisant la fonction projectRaster () dans le package raster en R. Un exemple reproductible complet est collé ci-dessous.

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

Ma question est la suivante: cet avertissement est-il un problème que je dois résoudre si je travaille sur des données terrestres? En d'autres termes, les données sont "perdues". Ce serait un problème majeur pour moi s'il l'est. Si tel est le cas, savez-vous s'il existe un moyen de résoudre ce problème?

J'ai cherché une solution à ce problème en ligne et j'en ai trouvé une mention ici , ici et ici, mais je pense qu'aucun ne propose une réponse appropriée à ce problème.

2. Chargez la bibliothèque raster

  library(raster)

3. Créez d'abord une carte (longue lat) du monde, avec une constante dans chaque cellule de la grille

Je mets une constante dans chaque cellule de la grille pour voir si je peux diagnostiquer le problème de l'endroit où l'avertissement est susceptible d'affecter les données, le cas échéant.

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

4. Si vous reprojetez la carte (lat longue) sur une grille de surface égale, elle produit le message d'avertissement

   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

Je pense que ce message dit essentiellement que la re-projection a échoué pour certaines des cellules de la grille.

5. Mais si vous tracez les deux cartes, il ne semble pas que cet avertissement pose des problèmes pour les données

Autrement dit, vous ne voyez pas de lacunes blanches dans les données tracées. Je suppose que les cellules perdues sont des cellules non terrestres, au sommet et au bord de l'étendue du monde. Des idées?

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)

entrez la description de l'image ici

Mike
la source

Réponses:

8

Réponse courte: ça va, et vous pouvez remercier d' rasteravoir procédé à l'achèvement au lieu d'un échec, et de vous avoir informé que certaines données ont été perdues.

Longue réponse:

cela dépendra de la projection, et dans ce cas c'est probablement juste sur les "bords". ce qu'est l'arête et comment elle se manifeste pour une instance donnée d'une famille de projection donnée est la partie "dépend".

Vous pouvez voir que ce ne sont pas les points centraux des cellules qui sont perdus:

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

Mais ce sont probablement les coins et peut-être les bords de chaque cellule. Cela montre peut-être que les projets raster sont basés sur l'étendue des cellules, et pas seulement sur leurs points centraux.

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

J'admets que je m'attendais à ce que ce soit d'où viennent les valeurs manquantes, alors peut projectRaster- être s'étend-il un peu plus au nord et au sud? Définissez-y des valeurs pour la latitude en dehors de la plage -90/90 et vous commencez à recevoir l'avertissement. Je ferai un suivi si j'ai l'occasion d'explorer davantage.

Enfin, vous devriez probablement utiliser un ellipsoïde explicite ou un paramètre de référence, c'est-à-dire "+ proj = eck4 + ellps = WGS84".

mdsumner
la source
Merci! Je viens de réinitialiser l'étendue du raster d'origine à moins / plus de +/- 70 degrés: rastertest.longlat <- raster (ncol = 4320, nrow = 2160, xmn = -180, xmx = 180, ymn = -70 , ymx = 70) puis reprojeté en utilisant: rastertest.eck4 <-projectRaster (rastertest.longlat, res = c (100000,100000), crs = "+ proj = eck4, method =" ngb ", over = T). l'avertissement disparaît complètement, donc je pense que toutes les valeurs manquantes doivent se produire +/- ~ 70 degrés.
Mike
@mike Je suis désolé mais ce n'est pas une solution acceptable, perdre des données et des informations n'est pas une solution. Qgis et Arcgis peuvent effectuer cette opération sans problème, il doit y avoir un problème avec projectRaster. J'espère que le développeur pourra jouer.
Herman Toothrot
(-‸ლ) Je vous invite à montrer exactement ce que font QGIS et Arc dans ce cas, spéculer n'est pas vraiment utile. Nous avons accès à exactement ce qui se passe dans R ici, il n'y a pas de problème réel.
mdsumner
5

Ce n'est pas une réponse complète à ma question initiale, en termes de tous les détails, mais cela fournit au lecteur intéressé un endroit où aller.

En général, les données sont perdues et déformées lors de la reprojection de rasters de projections longues à surfaces égales. Cela peut être problématique pour l'analyse globale. Si vous pouvez obtenir vos données au format vectoriel, il est préférable de reprojeter en utilisant ce format à la place.

Voici une référence sur le problème général. Et voici un autre sur la quantification de la perte. J'espère que cela pourra aider.

Mike
la source