У меня есть набор данных по сетке километров в континентальной части США. Столбцы «широта», «долгота» и «наблюдение», например:
"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
или, как фрейм данных 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))
(полный набор данных можно скачать как CSV здесь )
Данные выводятся из модели культуры (предназначенной для включения) в сетку 30 км х 30 км (из Miguez et al 2012 ).
Как я могу преобразовать их в растровый файл с метаданными, связанными с ГИС, такими как проекция карты?
В идеале файл должен быть текстовым (ASCII?), Потому что я хотел бы, чтобы он был независимым от платформы и программного обеспечения.
proj +proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +lon_0=-145.5 +lat_0=1.0
. Я не уверен, что вы запрашиваете в отношении другого CSV-файла - как он будет отличаться от того, на который ссылается вопрос, или что будет получено в результате принятого ответа?Ответы:
Требуется несколько шагов:
Вы говорите, что это обычная 1-километровая сетка, но это означает, что широта не является регулярной. Сначала вам нужно преобразовать его в регулярную систему координат сетки, чтобы значения X и Y регулярно располагались.
а. Прочитайте его в R как фрейм данных со столбцами x, y, yield.
б. Преобразуйте фрейм данных в SpatialPointsDataFrame, используя пакет sp и что-то вроде:
с. Преобразуйте в свою обычную систему km, сначала сообщив ей, что это за CRS, а затем spTransform до места назначения.
д. Скажите R, что это сетка:
В этот момент вы получите ошибку, если ваши координаты не лежат на хорошей регулярной сетке.
Теперь используйте растровый пакет для преобразования в растр и установите его CRS:
Теперь посмотрим:
Теперь запишите его как файл geoTIFF, используя растровый пакет:
Этот geoTIFF должен быть доступен для чтения во всех основных пакетах ГИС. Очевидным недостающим фрагментом здесь является строка proj4 для преобразования: это, вероятно, будет какая-то система ссылок UTM. Трудно сказать без каких-либо дополнительных данных ...
источник
r = Import[ "https://ebi-forecast.igb.illinois.edu/bety/miscanthusyield.csv", "Data"];
. Вы можете получить быстрый график точек затем черезdata = Rest[r]; ListPlot[data[[;; , {3, 2}]]]
(илиListPointPlot3D[data[[;; , {3, 2, 4}]]]
). Для репроекций, начните с помощиGeoGridPosition
, затем сделайте некоторые умные предположения и перекрестные ссылки, чтобы выяснить, что происходит :-). Кстати, объяснение @ Spacedman действительно актуально: метрическое искажение от 25 до 49 градусов равно cos (25) / cos (49) = 1,38; это существенно.Поскольку вопрос был последний ответил, гораздо проще решение существует с помощью растра ПАКЕТ в
rasterFromXYZ
функцию , которая инкапсулирует все шаги , необходимые (включая спецификацию строки CRS).источник
/tmp/
ГБ, когда мне не хватило места на диске.merge()
результаты вместе.