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

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.

— 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:

Wymagane kilka kroków:

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.
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")
```
Teraz spójrz:
```
plot(r)
```
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

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

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

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

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

... 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