Jak przekonwertować dane w postaci wartości lat, lon, value na plik rastrowy za pomocą R?

40

Mam zestaw danych wartości ponad kilometrową siatką w kontynentalnych Stanach Zjednoczonych. Kolumny to „szerokość geograficzna”, „długość geograficzna” i „obserwacja”, np .:

"lat"    "lon"     "yield"
 25.567  -120.347  3.6 
 25.832  -120.400  2.6
 26.097  -120.454  3.4
 26.363  -120.508  3.1
 26.630  -120.562  4.4

lub, jako ramka danych R:

mydata <- structure(list(lat = c(25.567, 25.832, 26.097, 26.363, 26.63), 
lon = c(-120.347, -120.4, -120.454, -120.508, -120.562), 
yield = c(3.6, 2.6, 3.4, 3.1, 4.4)), .Names = c("lat", 
"lon", "yield"), class = "data.frame", row.names = c(NA, -5L))

(pełny zestaw danych można pobrać tutaj jako csv )

Dane są wyprowadzane z modelu uprawy (przeznaczonego do włączenia) na siatce 30 km x 30 km ( Miguez i in. 2012 ).

wprowadź opis zdjęcia tutaj

Jak przekonwertować je na plik rastrowy za pomocą metadanych związanych z GIS, takich jak odwzorowanie mapy?

Idealnie byłby to plik tekstowy (ASCII?), Ponieważ chciałbym, aby był niezależny od platformy i oprogramowania.

raster convert r Abe
źródło
Jako CSV, to już to jest „plik tekstowy” w ASCII. Ponadto, ponieważ w ogóle nie korzysta z projekcji, do dodania może być niewiele istotnych metadanych (głównie dane odniesienia). Czy mógłbyś bardziej szczegółowo określić, jakiego rodzaju wyników szukasz i co zamierzasz z tym zrobić?
whuber
Chciałbym ułatwić użytkownikom korzystanie z danych za pomocą różnych programów do mapowania (ArcGIS, Google Maps, Grass, R itp.), Aby ułatwić ich ponowne użycie, np. Nie wymagając dodatkowych kroków konwersji. Na podstawie strony Wikipedii dotyczącej formatów plików GIS wnioskuję: 1) plik „rastrowy” powinien mieć nazwy geograficzne o szerokości i nazwach kolumn o długości geograficznej, np. Obraz, oraz że 2) metadane powinny zawierać informacje geograficzne (położenie rogu, obszar objęty według danych).
Abe
Jest to jedna z najlepszych referencji, na jakie natknąłem się na R i GIS. Dziękuję Ci bardzo! Czy możesz podać inny plik CSV z łatą Lat i Long z prawidłowym proj4string? Naprawdę to docenię.
@Nandini Nie wiem, co jest poprawne proj4string, podejrzewam Lamberta:proj +proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +lon_0=-145.5 +lat_0=1.0 . Nie jestem pewien, o co prosisz w odniesieniu do innego pliku csv - czym różni się on od tego, do którego link znajduje się w pytaniu lub czy zostałby uzyskany na podstawie zaakceptowanej odpowiedzi?
Abe
dla mnie nie działa! Nie wiem, co umieścić na współrzędnych „x” i „y” na „(pts) = ~ x + y”

Odpowiedzi:

44

Wymagane kilka kroków:

  1. Mówisz, że jest to regularna siatka o długości 1 km, ale to oznacza, że ​​długość lat nie jest regularna. Najpierw musisz przekształcić go w regularny układ współrzędnych siatki, aby wartości X i Y były regularnie rozmieszczone.

    za. Wczytaj to do R jako ramkę danych, z kolumnami x, y, wydajnością.

    pts = read.table("file.csv",......)

    b. Konwertuj ramkę danych na SpatialPointsDataFrame przy użyciu pakietu sp i czegoś takiego:

    library(sp)
library(rgdal)
coordinates(pts)=~x+y

    do. Konwertuj na swój zwykły system km, najpierw informując go, co to jest CRS, a następnie spTransformuj do miejsca docelowego.

    proj4string(pts)=CRS("+init=epsg:4326") # set it to lat-long
pts = spTransform(pts,CRS("insert your proj4 string here"))

    re. Powiedz R, że jest to gridowane:

    gridded(pts) = TRUE

    W tym momencie pojawi się błąd, jeśli twoje współrzędne nie leżą na ładnej regularnej siatce.

  2. Teraz użyj pakietu rastrowego, aby przekonwertować go na raster i ustawić jego CRS:

    r = raster(pts)
projection(r) = CRS("insert your proj4 string here")

  3. Teraz spójrz:

    plot(r)

  4. Teraz napisz go jako plik geoTIFF, używając pakietu rastrowego:

    writeRaster(r,"pts.tif")

Ten geoTIFF powinien być czytelny we wszystkich głównych pakietach GIS. Oczywistym brakującym elementem jest ciąg proj4, na który należy przekonwertować: prawdopodobnie będzie to jakiś system odniesienia UTM. Trudno powiedzieć bez dodatkowych danych ...

Spacedman
źródło
+1 Dzięki za ułożenie przepływu pracy. Pamiętaj, że dane są dostępne pod linkiem podanym w pytaniu: spójrz. Odkryjesz niestety, że niektóre z twoich założeń na ich temat są błędne. (W szczególności, polowałem na każdej dokumentacji dotyczącej projekcji użyte do wygenerowania siatki, ale nic nie znaleźli I to jest dziwne, projekcja, jak widać przez wykreślenie pkt.).
whuber
Jest bardzo bliski bycia systemem UTM, ale żaden z tych, których próbowałem, nie jest wystarczająco blisko zwykłej siatki, aby R mógł je obrysować. Jestem w połowie kuszony, aby przejrzeć całą bazę danych epsg R. ....
Spacedman,
To byłby prawdziwy tour de force, gdybyś mógł odkryć projekcję w ten sposób! Kluczem jest znalezienie skutecznego i wydajnego kryterium, aby ustalić, kiedy te ponad 7000 punktów jest wystarczająco blisko, aby leżeć na regularnej siatce (ponieważ jest możliwe, że mogą nie tworzyć idealnej siatki w żadnej standardowej projekcji). Aby szybko przeszukać bazę danych, powinno wystarczyć porównanie niewielkiej liczby odległości, takich jak odległość wschód-zachód w północnej części siatki do odległości wschód-zachód w części południowej. To powinno szybko wyeliminować ogromną większość kandydatów.
whuber
3
Przejrzałem wszystkie (domyślne) projekcje obsługiwane przez Mathematica 8. Znalazłem projekcję, w której punkty naprawdę wydają się spadać na siatkę: Alaska State Plane (1983) Zone 10! Jest to projekcja stożkowa Lamberta. Myślę, że to EPSG 26940 . Jeśli zmodyfikujesz to, aby wyśrodkować go w przybliżeniu na długości -106, punkty tworzą całkiem dobrą siatkę.
whuber
1
Abe, masz na myśli przeczytać stronę internetową? Tak było r = Import[ "https://ebi-forecast.igb.illinois.edu/bety/miscanthusyield.csv", "Data"];. Następnie możesz uzyskać szybki wykres punktów za pomocą data = Rest[r]; ListPlot[data[[;; , {3, 2}]]](lub ListPointPlot3D[data[[;; , {3, 2, 4}]]]). W przypadku ponownych rzutów zacznij od pomocy GeoGridPosition, a następnie zrób inteligentne domysły i odsyłacze, aby dowiedzieć się, co się dzieje :-). BTW, wyjaśnienie @ Spacedmana jest naprawdę istotne: zniekształcenie metryczne od 25 do 49 stopni jest równe cos (25) / cos (49) = 1,38; to znaczące.
whuber
29

Od czasu ostatniej odpowiedzi na pytanie istnieje znacznie łatwiejsze rozwiązanie, używając funkcji pakietu rastrowego, rasterFromXYZktóra zawiera wszystkie niezbędne kroki (w tym specyfikację ciągu CRS). 

library(raster)
rasterFromXYZ(mydata)
Lucas Fortini
źródło
1
Przepraszam niestrudzonego @Spacedmana, który często mi pomagał, ale myślę, że ta odpowiedź zasługuje na odziedziczenie wesołego zielonego kleszcza.
geotheory
@geotheory Wybrałbym tę odpowiedź, jest to świetna funkcja, ale wydaje się, że jest bardzo powolna w zestawie danych, którego używam (połączonym w op)
Abe
1
... w rzeczywistości zadławił się, ponieważ wziął mój plik ~ 400 KB i utworzył plik o rozmiarze ~ 19 GB, /tmp/kiedy zabrakło mi miejsca na dysku.
Abe
Prawdopodobnie jest gdzieś tam proces n-kwadrat. Możesz pogrupować dane punktów za pomocą szerokiej siatki, zrasteryzować każdą grupę osobno, a następnie merge()wyniki razem.
geotheory
Z całym szacunkiem, ale ta odpowiedź jest znacznie lepsza niż odpowiedź Spacedmana.
Ghost