Измените crs растра, чтобы он соответствовал crs простого точечного объекта
Допустим, у меня есть растр MODIS LAI и простой объект (sf) (точки) с двумя разными crs. Мне нужно преобразовать crs растра, чтобы он соответствовал crs данных sf (и наоборот, тоже можно, но я бы предпочел избавиться от синусоидальной системы MODIS). Различные попытки не привели меня к решению, т. е. иметь растр и sf в одной CRS для дальнейшей визуализации, обработки и так далее. Любая помощь о том, как я могу это сделать, пожалуйста? Ниже приведены несколько вариантов, которые я попробовал:
library(terra)
library(sf)
lai <- terra::rast("lai.tif")
terra::crs(lai) #check crs of lai
plots <- readRDS("plots.rds")
sf::st_crs(plots) #check crs of sf
#Trial 1
plots_tr <- sf::as_Spatial(plots) #transform the sf and get the crs
lai_pr1 <- terra::project(lai, "+proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") #project the raster using the crs of the transformed sf
terra::plot(lai_pr1) #plot raster
plot(plots_tr, add=T, pch=19, col = "black", cex=0.1) #overlay points
#Trial 2
lai_pr2 <- lai
terra::crs(lai_pr2) <- "+proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" # replace lai crs
terra::plot(lai_pr2) #plot raster
plot(plots_tr, add=T, pch=19, col = "black", cex=0.1) #overlay points
#Trial 3
lai_pr3 <- terra::project(lai, "epsg:25831") #project the raster using the EPSG code of the sf data (i.e., ETRS89 / UTM zone 31N)
terra::plot(lai_pr3) #plot raster
plot(plots_tr, add=T, pch=19, col = "black", cex=0.1) #overlay points
Ссылка на данные для запуска этого примера:
лай <- ссылка удалена
сюжеты <- ссылка удалена
Может быть, проблема в том, что crs sf — это WKT, а crs растра — нет? Кроме того, использование sf без sf::as_Spatial() не изменило результат. Я также попробовал изменить crs научной фантастики с помощью st_transform(). Я прочитал множество вопросов, уже заданных по этому поводу, но не смог найти решения. Maybe spTransform()?
Надеюсь, это не слишком глупо. Большое спасибо.
Кажется, это та же проблема, что и на днях, когда у вас возникли проблемы, связанные с crs, с данными MODIS. Происходит сдвиг между растровыми и векторными данными. ГЕ снова? terra::crs(lai) <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"
Ответы 1
После исправления системы координат вашего файла MODIS это довольно просто, используя sf::st_transform() или terra::project(), в зависимости от того, какие crs вы хотите сохранить. Для анализа я бы продолжил перепроецирование векторных данных и оставил бы растр как есть, но не стесняйтесь избавиться от проекции MODIS, если хотите (только не забывайте о повторной выборке). Также совершенно не обязательно использовать {sp}.
library(terra)
#> terra 1.7.71
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
lai <- terra::rast("lai.tif")
plots <- readRDS("plots.rds")
# fix MODIS crs
modis_crs <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"
terra::crs(lai) <- modis_crs
# option 1: reproject vector data
plots_reproj <- sf::st_transform(plots, modis_crs)
terra::plot(lai, main = "MODIS Sinusoidal")
plot(plots_reproj, add = TRUE, pch = 19, col = "black", cex = 0.1)
# option 2: reproject raster data
lai_reproj <- terra::project(lai, "epsg:25831")
terra::plot(lai_reproj, main = "EPSG:25831")
plot(plots, add = TRUE, pch = 19, col = "black", cex = 0.1)
Created on 2024-04-27 with reprex v2.1.0
Да, мой вопрос был глупым :). Дело в том, что это вычисление происходит последовательно после проблемы , которая у меня была раньше . Когда эта проблема была решена, я подумал, что crs MODIS в порядке, поэтому я не уверен, почему crs MODIS снова ошибся (но я мог бы проверить). Имеет смысл? CRS теряется каждый раз, когда я выполняю вычисления? В любом случае, большое спасибо за вашу дальнейшую помощь. Вы также получаете предупреждающее сообщение после запуска sf::st_transform(plots, modis_crs)? Спасибо!
Пространственные атрибуты сохраняются и не должны теряться во время вычислений и т. д. Я предполагаю, что вы, возможно, исправили это для какого-то объекта SpatRaster, но все становится запутанным, и вы, вероятно, переходите к работе с другим. Или, может быть, просто перезаписать предыдущий (фиксированный) объект объектом, импортированным из (нефиксированных) файлов еще раз? Лучше проверьте свой код, такого быть не должно. Кроме того, мне не вернулось никаких предупреждений. Что это говорит?
Да, я еще раз проверю свой код, спасибо, так как согласен с вами, что этого не должно происходить. В соответствии с предупреждением я получаю следующее: > Предупреждающее сообщение: > В CPL_transform(x, crs, aoi, конвейер, обратный, желаемая_точность,: > Ошибка GDAL 1: PROJ: proj_as_wkt: DatumEnsemble можно экспортировать только в > WKT2:2019
@Фрэн, иногда гуглишь «неясные последние части ошибки», здесь «DatumEnsemble можно экспортировать только в >WKT2:2019», это может привести к подсказкам StackOverflow приводит к r-spatial sf и ты можешь разрабатывать диагноз: я использую Windows или Linux (Mac), какой у меня дистрибутив, не слишком ли старый мой PROJ GDAL GEOS, особенности, исправления ошибок и тому подобное. Но помогают именно последние, странные слова, которые по большей части уникальны и приводят к интересному чтению.
?terra::project, вы увидите, что рекомендуется, чтобы буква «y» в project(x, y также была Spatrast, также была SpatRast, поэтому растрируйте свой «sf», а затем используйте project. Посмотрим немного ближе, но это общую идею смотрите Ответ Космонавта.