R: Jak uzyskać szerokości i długości geograficzne z warstwy RasterLayer?

14

Jestem absolutnym początkującym użytkownikiem danych geograficznych, więc proszę wybacz mi, jeśli pytanie nie jest właściwe.

Pobrałem dane z NCDC NARR i udało mi się załadować do R przy użyciu rasterpakietu. Chciałbym uzyskać listę z szerokością, długością i wartością. Rozumiem, że rasterToPoints()powinienem robić dokładnie to, co chcę, jednak moje wartości szerokości i długości geograficznej wyglądają dziwnie:

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

Przypuszczam, że powinienem coś zrobić z projekcją, którą jest obecnie Lambert Conformal Conic (LCC). Oto dalsze informacje o rastrze.

> 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

Co mam zrobić, aby uzyskać rzeczywiste szerokości i długości geograficznej w USA?

janosdivenyi
źródło

Odpowiedzi:

14

musisz zmienić projekt rastra na rzut geograficzny (stopnie dziesiętne) za pomocą „projectRaster” lub „spTransform”. Zobacz także definicje sp CRS, które określają pożądany ciąg projekcji. Przykład w pomocy dla „projectRaster” jest dość jasny, jak to zrobić.

Jeśli zmusisz swoje dane rastrowe do obiektu SpatialPointsDataFrame, wówczas użyjesz „spTransform” i wyciągniesz współrzędne ze szczeliny @coordinates i dodasz je do data.frame w szczelinie @data. Oto przykład tego, jak by to wyglądało.

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)

Powinienem zauważyć, że konwersja rastrów na klasę obiektów wektorowych nie jest dobrą praktyką i neguje zalety pakietu rastrowego zapewniającego bezpieczne przetwarzanie pamięci. Często rozsądnie jest myśleć o swoim problemie i oceniać, czy podchodzisz do niego poprawnie. Jeśli PO podał kontekst, dlaczego potrzebują współrzędnych [x, y] dla każdej komórki, społeczność forum mogłaby być w stanie zapewnić alternatywne obliczenia, które utrzymałyby problem w środowisku rastrowym.

Jeffrey Evans
źródło
1
Jednym ze sposobów na zachowanie ostrożności (unikanie konwersji danych) do serca jest odznaczenie pierwotnego rastra (być może na bardzo grubej siatce), utworzenie dwóch siatek wartości szerokości i długości geograficznej obejmujących zakres nie projekcji i rzutowanie ich z powrotem na zasięg oryginalnej siatki. Nie są tworzone klasy wektorowe: jest to całkowicie zestaw operacji rastrowych.
whuber
4

Uzyskaj współrzędne centrów komórek i utwórz obiekt przestrzenny:

spts <- rasterToPoints(r, spatial = TRUE)

Przekształć punkty w pożądany cel:

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

Wartości są już skopiowane jako kolumny w tym SpatialPointsDataFrame.

print(llpts)

Teraz, aby zakończyć, zdobądź data.frame:

x <- as.data.frame(llpts)

Istnieje ogólna implementacja tego w pakiecie SGAT, patrz funkcja lonlatFromCelltutaj:

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

mdsumner
źródło
Próbowałem tego, ale otrzymałem następujący komunikat o błędzie: > 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
Właściwie nie musisz przenosić atrybutów, usunę to.
mdsumner,
Oprócz porady dotyczącej pakietu SGAT, czy nie jest to dokładnie ta sama odpowiedź / przykład, który podałem? Współrzędne nie są propagowane do data.frame w gnieździe danych, tylko wartości z rastra. Współrzędne są w rzeczywistości przechowywane w szczelinie współrzędnych i należy je dodać do data.frame.
Jeffrey Evans,
Dzięki, dodałem krok as.data.frame. Myślę, że dodawanie współrzędnych jako atrybutów jest okropną wskazówką - szczególnie poprzez mungowanie ze szczeliną - ponieważ współrzędne obiektu mogą się zmieniać. Jeśli chcesz surową ramkę danych, po prostu ją zrób. Nie dbam o to, gdzie są informacje, może po prostu edytuj swoje i możemy zapełnić tę odpowiedź.
mdsumner,
OP konkretnie chciał współrzędnych i myślę, że zapisywanie w osobnej ramce danych jest zbędne. Zwykle nie lubię dodawać współrzędnych głównie do szczeliny danych, ponieważ jest ona zbędna w stosunku do szczeliny współrzędnych. Poza tym dodawanie informacji do szczeliny danych nie jest „okropną radą”. Co jeśli chcesz mieć dwa układy współrzędnych. Możesz dodać długość / długość do szczeliny danych i ustawić obiekt w zupełnie innej projekcji. Ponadto, jeśli chcesz po prostu wyeksportować płaski plik, a nie format GIS jako taki, możesz dodać współrzędne do data.frame i zapisać jako plik csv.
Jeffrey Evans,
0

Wygląda na to, że masz tam rzutowane współrzędne (nie szerokość / długość geograficzna, czyli współrzędne GCS). Prawdopodobnie nie było dla ciebie jasne, że to był problem. Zobacz ten post. Przekształcanie układu współrzędnych geograficznych w R.

jbchurchill
źródło
Nie złapałem posta, do którego odwoływałeś się, zanim odpowiedziałem. Możesz oznaczyć to jako duplikat. Chociaż dodanie przymusu i przypisania współrzędnych SpatialPointsDataFrame sprawia, że ​​jest nieco inaczej. Twoja decyzja.
Jeffrey Evans,
Myślałem o oznaczeniu go jako takiego, ale pomyślałem, że jeśli inna osoba szuka podobnej odpowiedzi, nie wiedząc, że musi rzutować wartości, może się to okazać. Poza tym twoja odpowiedź oferuje inny sposób dotarcia do tego miejsca (pozytywny).
jbchurchill 10.04.15
Próbowałem spojrzeć na wymienione przez ciebie źródła. W celu uzyskania standardowych szerokości / długości geograficznych wydałem lonlat_r <- projectRaster(r, crs="+init=epsg:4326"). Jednak zakres nowego rastra jest -181.3232, 181.4938, -1.590457, 87.76154 (xmin, xmax, ymin, ymax)daleki od tego, czego oczekiwałbym od USA (który powinien wynosić od 30 do 70 i od -60 do -160). Powinienem coś źle zrozumieć.
janosdivenyi