Значения RGB упали после проецирования SpatRaster
Как указано в заголовке, после запуска project() в файле .tif со слоями RGB значения отбрасываются, и при построении графика возвращаются только NA. Ближайшее к решению, которое я нашел, это здесь и здесь, но в обоих случаях значения сохранились после project() (хотя и в измененной форме).
Я пробовал различные project() варианты, в том числе проецирование отдельных слоев, но безрезультатно. Я подозреваю, что требуется другой тип класса? Однако я не смог найти ничего в документации или каких-либо окончательных отчетов об ошибках GitHub.
Репекс:
library(basemaps)
library(terra)
library(sf)
# Create sf for SpatRaster extent
df <- data.frame(lon = c(170, 172),
lat = c(-43, -45))
sf_poly <- df %>%
st_as_sf(coords = c("lon", "lat"),
crs = 4326) %>%
st_bbox() %>%
st_as_sfc() %>%
st_as_sf() %>%
st_transform(3857)
r <- basemap_raster(ext = st_bbox(sf_poly)) %>%
as("SpatRaster")
r
# class : SpatRaster
# dimensions : 1012, 728, 3 (nrow, ncol, nlyr)
# resolution : 305.7481, 305.7481 (x, y)
# extent : 18924279, 19146864, -5621485, -5312068 (xmin, xmax, ymin, ymax)
# coord. ref. : WGS 84 / Pseudo-Mercator
# source : basemap_20240418094850.270686.tif
# colors RGB : 1, 2, 3
# names : red, green, blue
# min values : 167, 161, 162
# max values : 255, 255, 255
Метод 1: Полный проект SpatRaster
r4326 <- project(r, "EPSG:4326")
r4326 <- project(r, "+proj=longlat +datum=WGS84 +no_defs +type=crs")
r4326 <- project(r, crs("EPSG:4326"))
r4326 <- project(r, "EPSG:4326", method = "bilinear")
r4326 <- project(r, "EPSG:4326", method = "near")
r4326 <- project(r, "EPSG:4326", align = TRUE)
r4326 <- project(r, "EPSG:4326", mask = TRUE)
r4326
# class : SpatRaster
# dimensions : 881, 882, 3 (nrow, ncol, nlyr)
# resolution : 0.002268062, 0.002268062 (x, y)
# extent : 169.9997, 172.0001, -44.99879, -43.00063 (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326)
Метод 2: Проецирование отдельных слоев SpatRaster
r_list <- split(r, 1:nlyr(r)) %>% `names<-`(names(r))
r_list[[1]]
# class : SpatRaster
# dimensions : 1012, 728, 1 (nrow, ncol, nlyr)
# resolution : 305.7481, 305.7481 (x, y)
# extent : 18924279, 19146864, -5621485, -5312068 (xmin, xmax, ymin, ymax)
# coord. ref. : WGS 84 / Pseudo-Mercator
# source : basemap_20240418094850.270686.tif
# name : red
# min value : 167
# max value : 255
r1 <- r_list[[1]] %>% project("EPSG:4326")
r1
# class : SpatRaster
# dimensions : 881, 882, 1 (nrow, ncol, nlyr)
# resolution : 0.002268062, 0.002268062 (x, y)
# extent : 169.9997, 172.0001, -44.99879, -43.00063 (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326)
@dimfalk, воспроизводимый задает вопрос, особенность или ...? и, вероятно, Роберту будет интересно. Я следил за тем, как мои попытки обычно не завершаются с первой попытки, как указано выше, но, возможно, я не был настолько скрупулезным в подсчете «один», «два». Кажется, стоит отчитаться. Возможно, предполагаемый «writeRaster» после проекта и не полагающийся на объект, но я размышляю о нем.
Ответы 1
До сих пор не знаю, почему это не получается, но, по крайней мере, я мог бы сузить это до basemap_raster(ext = st_bbox(sf_poly)) %>% as("SpatRaster"). Может быть, некоторые временные файлы больше не доступны после некоторой обработки? Просто предполагаю.
Но если вместо этого вы используете basemap_terra(), все ваши версии project() будут предоставлять значения RGB без каких-либо NA.
library(basemaps)
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
# Create sf for SpatRaster extent
df <- data.frame(lon = c(170, 172),
lat = c(-43, -45))
sf_poly <- df %>%
st_as_sf(coords = c("lon", "lat"),
crs = 4326) %>%
st_bbox() %>%
st_as_sfc() %>%
st_as_sf() %>%
st_transform(3857)
r <- basemap_terra(ext = st_bbox(sf_poly))
#> Loading basemap 'voyager' from map service 'carto'...
r
#> class : SpatRaster
#> dimensions : 1012, 728, 3 (nrow, ncol, nlyr)
#> resolution : 305.7481, 305.7481 (x, y)
#> extent : 18924279, 19146864, -5621485, -5312068 (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
#> source : basemap_20240418151816.750852.tif
#> colors RGB : 1, 2, 3
#> names : red, green, blue
#> min values : 65, 92, 120
#> max values : 254, 248, 246
r4326 <- project(r, "EPSG:4326")
r4326
#> class : SpatRaster
#> dimensions : 881, 882, 3 (nrow, ncol, nlyr)
#> resolution : 0.002268062, 0.002268062 (x, y)
#> extent : 169.9997, 172.0001, -44.99879, -43.00063 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source(s) : memory
#> colors RGB : 1, 2, 3
#> names : red, green, blue
#> min values : 65, 92, 120
#> max values : 254, 248, 246
Совершенно пропустил library(basemap) выше. Интересное первое погружение здесь, так как as('SpatRaster') должен явно указать цель terra (или я бы так предположил).
@Chris: Это также не работает с terra::rast() вместо этого и без трубопровода, так что, может быть, что-то с объектом растрового кирпича, возвращаемым basemap_raster(), не так?
Я думаю, мы упростим задачу, исключив basemap из цепочки (поскольку до сих пор у меня его не было), и поэтому моя подписка/комментарий основывалась просто на terra и задавалась вопросом, почему, как и @l-tyrone, мои ценности в проекте, похоже, упали. , что, как я предполагал, было чем-то, что я делал неправильно, но предполагал, что значения будут применяться ко всей проецируемой цели. Тем не менее, это может быть просто процедурный вопрос, который нужно решить в ?, или это может быть функция базового кода. Вы заметили это, ваш отчет.
Черт, похоже, это не способ получить 1111, вам придется редактировать очередь больше, чтобы изучить нечетные числа.
basemap, в их readme подразумевается эквивалентность вызовов между basemap_raster() и basemap_terra(), но я думаю, что это все еще работает с terra.
@Chris: Чтобы сохранить это минимумом, мы могли бы также использовать system.file("ex/logo.tif", package = "terra") |> rast() в качестве растра RGB, но я не могу воспроизвести указанное поведение, используя эти данные... Вот почему я был все более уверен, что в конце концов речь идет о {basemaps}. Все еще интересно, по какой причине это не помогло вам.
Основываясь на вашем ответе, я нашел источник проблемы. Вы, @Chris, и больше всего я будем ругать себя: basemap_raster() требует загрузки пакета raster! Учитывая, что вам не удалось воспроизвести проблему, я считаю, что проблема может заключаться в привязках GDAL, GEOS и PROJ. Я почти уверен, что мне пришлось внести некоторые исправления в тот же день. В любом случае, проголосуйте за и примите указание пути.
@LTyrone: Потрясающе! :) Позвольте мне упомянуть этот крошечный факт в ответе - на случай, если кто-то не прочитает все комментарии здесь.
@dimfalk - надеялся, что ты это сделаешь :)
@LTyrone: Я как раз собирался это сделать, но почему-то у меня это все еще не получается, даже с library(raster) до basemaps::basemap_raster(). То же поведение - срабатывает один раз, но не дважды.
@dimfalk - думаю, пора вообще отказаться от basemap_raster(). У меня сейчас тоже со второй попытки не получается. Я отправлю отчет об ошибке.
Я собирался написать, что не могу воспроизвести это, потому что мой результат
project(r, "EPSG:4326")выглядел вполне нормально, но потом после первой успешной попытки это внезапно начало давать сбой. Перезапуск R позволил мне перепроецировать растр один раз, но никогда во второй раз без потери значений. Это дико. Ошибка?