R: Comment obtenir les latitudes et longitudes d'un RasterLayer?

14

Je suis un débutant absolu des données géographiques, alors s'il vous plaît, pardonnez-moi si la question n'est pas appropriée.

J'ai téléchargé des données depuis NCDC NARR et j'ai réussi à charger dans R à l'aide du rasterpackage. Je voudrais obtenir une liste avec latitude, longitude et valeur. Je comprends que cela rasterToPoints()devrait faire exactement ce que je veux, cependant, mes valeurs de latitude et de longitude semblent étranges:

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

Je suppose que je devrais faire quelque chose avec la projection qui est actuellement la conique conforme de Lambert (LCC). Voici d'autres informations sur le raster.

> 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

Que dois-je faire pour obtenir des valeurs réelles de latitude et de longitude américaines?

janosdivenyi
la source

Réponses:

14

vous devez réellement reprojeter le raster dans une projection géographique (degrés décimaux) à l'aide de "projectRaster" ou "spTransform". Regardez également les définitions de CRS sp qui spécifient la chaîne de projection souhaitée. L'exemple dans l'aide du "projectRaster" est assez clair sur la façon de procéder.

Si vous contraignez vos données raster dans un objet SpatialPointsDataFrame, vous utiliserez alors "spTransform" et tirerez les coordonnées de l'emplacement @coordinates et les ajouterez au data.frame dans l'emplacement @data. Voici un exemple de ce à quoi cela ressemblerait.

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)

Je dois noter que ce n'est pas une bonne pratique de convertir des rasters en une classe d'objets vectoriels et que cela annule les avantages du package raster offrant un traitement sûr de la mémoire. Il est souvent prudent de vraiment réfléchir à votre problème et d'évaluer si vous l'abordez correctement. Si l'OP avait fourni un contexte expliquant pourquoi ils avaient besoin des coordonnées [x, y] pour chaque cellule, la communauté du forum aurait pu être en mesure de fournir des alternatives de calcul qui permettraient de conserver le problème dans un environnement raster.

Jeffrey Evans
la source
1
Une façon de prendre à cœur votre prudence (concernant l’évitement de la conversion des données) consiste à dé-projeter le raster original (peut-être sur une grille très grossière), à ​​créer deux grilles de valeurs de latitude et de longitude couvrant l’étendue de la non-projection et à les projeter à nouveau dans le étendue de la grille d'origine. Aucune classe vectorielle n'est créée: il s'agit entièrement d'un ensemble d'opérations raster.
whuber
4

Obtenez les coordonnées des centres cellulaires et créez un objet spatial:

spts <- rasterToPoints(r, spatial = TRUE)

Transformez les points en cible souhaitée:

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

Les valeurs sont déjà copiées sous forme de colonnes sur ce SpatialPointsDataFrame.

print(llpts)

Maintenant, pour terminer, obtenez un data.frame:

x <- as.data.frame(llpts)

Il y a une implémentation générale de cela dans le package SGAT, voir la fonction lonlatFromCellici:

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

mdsumner
la source
J'ai essayé mais j'ai reçu le message d'erreur suivant: > 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
En fait, vous n'avez pas besoin de transférer les attributs, je vais les supprimer.
mdsumner
À part les conseils sur le package SGAT, n'est-ce pas exactement la même réponse / exemple que j'ai fourni? Les coordonnées ne sont pas propagées au data.frame dans la fente de données, seulement les valeurs du raster. Les coordonnées sont, en fait, conservées dans l'emplacement des coordonnées et doivent être ajoutées au data.frame.
Jeffrey Evans
Merci, j'ai ajouté l'étape as.data.frame. Je pense que c'est un mauvais conseil d'ajouter les coordonnées en tant qu'attributs - en particulier en les accouplant avec une fente - car les coordonnées de l'objet peuvent changer. Si vous voulez un data.frame brut, faites-en un. Peu m'importe où se trouvent les informations, modifiez peut-être les vôtres et nous pouvons zapper cette réponse.
mdsumner
L'OP voulait spécifiquement des coordonnées et je pense qu'il est redondant d'enregistrer sur un data.frame séparé. Normalement, je n'aime pas ajouter de coordonnées à l'emplacement de données principalement, car il est redondant avec l'emplacement de coordonnées. En dehors de cela, ce n'est pas un "mauvais conseil" d'ajouter des informations à la fente de données. Et si vous voulez avoir deux systèmes de coordonnées. Vous pouvez ajouter lat / long à la fente de données et avoir l'objet dans une projection entièrement différente. De plus, si vous souhaitez simplement exporter un fichier plat, et non un format SIG en soi, vous pouvez ajouter des coordonnées au data.frame et l'enregistrer au format csv.
Jeffrey Evans
0

Il semble que vous ayez des coordonnées projetées là-bas (pas Latitude / Longitude alias coordonnées GCS). Ce n'était probablement pas clair pour vous que c'était le problème. Voir cet article. Conversion du système de coordonnées géographiques dans R

jbchurchill
la source
Je n'ai pas attrapé le poste auquel vous faites référence avant de répondre. Vous voudrez peut-être le signaler comme doublon. Bien que l'ajout de la coercition et de l'affectation des coordonnées SpatialPointsDataFrame fasse un peu différent. Ton appel.
Jeffrey Evans
J'ai pensé à le marquer comme tel, mais j'ai pensé que si une autre personne cherche une réponse similaire sans savoir qu'elle doit projeter les valeurs, cela pourrait se révéler pour elle. De plus, votre réponse offre un moyen différent d'y arriver (vote positif).
jbchurchill
J'ai essayé de regarder les sources que vous avez énumérées. Afin d'obtenir des latitudes / longitudes standard, j'ai émis lonlat_r <- projectRaster(r, crs="+init=epsg:4326"). Cependant, l'étendue du nouveau raster est -181.3232, 181.4938, -1.590457, 87.76154 (xmin, xmax, ymin, ymax)loin de ce que j'attendrais des États-Unis (qui devrait se situer entre 30 et 70 et -60 à -160). J'aurais dû mal comprendre quelque chose.
janosdivenyi