Почему terra выдает разные диапазоны значений для одного и того же растра?
Я использую terra для чтения и последующего сохранения различных глобальных наборов растровых данных. Сегодня я заметил совершенно разные диапазоны значений в зависимости от того, как я смотрю на данные и что с ними делаю. Я сбит с толку и нервничаю, пытаясь провести какой-либо анализ, пока не смогу объяснить, почему я вижу то, что вижу. Я был бы очень признателен, если бы кто-нибудь из вас рассказал мне, почему диапазоны значений, которые я вижу, настолько различаются.
В качестве примера можно использовать глобальную плотность дорог. Я импортирую растр GRIP Total Density, все типы вместе ascii, используя terra следующим образом:
roadDir <- c("downloads/GRIP4_density_total/")
tmp <- paste0(roadDir, "grip4_total_dens_m_km2.asc")
roadDensity <- terra::rast(tmp)
Если я затем посмотрю на roadDensity, то увижу, что диапазон значений составляет от 0 до 99445 (м на квадратный километр):
> roadDensity
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source : grip4_total_dens_m_km2.asc
name : grip4_total_dens_m_km2
min value : 0
max value : 99445
Если я импортирую исходные данные в QGIS, я получаю тот же размер ячейки, размеры и диапазон значений, как показано в приведенном выше фрагменте кода, поэтому я принимаю эти значения как истинные значения в данных. Все идет нормально.
Проблемы возникают, когда: а) я просматриваю сводные данные или б) записываю полученный растр в файл .tif для использования в дальнейшей обработке.
а) Резюме
terra::global(roadDensity, c( "min", "max"), na.rm=TRUE)
min max
grip4_total_dens_m_km2 0 330306
Я не понимаю, почему теперь я получаю совсем другое максимальное значение (330306).
б) Напишите в .tif
x <- writeRaster(roadDensity,"test.tif", datatype = "INT4U", overwrite=TRUE)
x
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : test.tif
name : grip4_total_dens_m_km2
min value : 0
max value : 330306
Таким образом, global и writeRaster terra дают одинаковые результаты.
Я понимаю, что путем изменения значения типа данных на одну из его возможностей («INT1U», «INT2U», «INT4U») writeRaster меняет диапазон значений из-за максимально возможных значений, которые может содержать конкретная структура данных. Чего я не понимаю, так это того, как «INT4U» может так сильно увеличить диапазон значений по сравнению с оригиналом без какого-либо изменения размера ячейки. Кроме того, значение INT*U не должно изменять сводные значения, создаваемые global.
Я уверен, что это очень простая вещь, и у меня настал момент «ага», но я был бы признателен за любые рекомендации, которые помогут моему пониманию и, следовательно, аналитической уверенности.
Спасибо
Когда я запускаю roadDensity <- rast("grip4_total_dens_m_km2.asc"), затем global(roadDensity, c( "min", "max"), na.rm=TRUE), min == 0 и max == 330306. Проверьте свою среду R на наличие случайного SpatRaster или запустите новый сеанс R и повторите попытку.
Чтобы уточнить, после запуска roadDensity <- rast("grip4_total_dens_m_km2.asc"), а затем roadDensity не отображаются значения минимального/максимального значения по умолчанию. Мне нужно запустить roadDensity[] <- values(roadDensity), чтобы они отобразились в метаданных. Опять же, min == 0, max == 330306. Так что да, попробуйте запустить процесс снова с самого начала.
Какой код вы используете, чтобы «посмотреть на roadDensity»?
Спасибо, ДиДжей. Мне известно о выборке, но у меня сложилось впечатление, что global использует все данные, если вы не используете весовую функцию. Рад, что исправился.
Спасибо, Л Тайрон. Запустил новый сеанс R и тот же результат. В моей системе (ubuntu 22.04) мне просто нужно «запустить» объект roadDensity, чтобы увидеть представленные мной результаты.
Что будет, если запустить roadDensity[] <- values(roadDensity), а затем roadDensity? Единственное, что я могу предложить, если здесь не будет ответа, — это поднять вопрос на GitHub.
@L Tyrone «Что произойдет, если вы запустите roadDensity[] <-values(roadDensity), а затем roadDensity?» Получил тот же результат. Спасибо за попытку.
Ответы 1
Я этого не вижу. Я предполагаю, что QGIS создал дополнительный файл с приблизительными экстремальными значениями.
Получить данные в чистую папку
dir.create("test")
setwd("test")
url <- "https://dataportaal.pbl.nl/downloads/GRIP4/GRIP4_density_total.zip"
zip <- basename(url)
if (!file.exists(zip)) {
download.file(url, zip, mode = "wb")
unzip(zip)
}
мин и макс не сообщаются для файла ascii (поскольку они не хранятся в формате файла)
library(terra)
#terra 1.7.78
r <- rast("grip4_total_dens_m_km2.asc")
r
#class : SpatRaster
#dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
#resolution : 0.08333333, 0.08333333 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
#source : grip4_total_dens_m_km2.asc
#name : grip4_total_dens_m_km2
Создайте новый SpatRaster, чтобы показать минимальные и максимальные значения:
r * 1
#class : SpatRaster
#dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
#resolution : 0.08333333, 0.08333333 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
#source(s) : memory
#varname : grip4_total_dens_m_km2
#name : grip4_total_dens_m_km2
#min value : 0
#max value : 330306
Или вычислите их следующим образом:
minmax(r, compute=TRUE)
# grip4_total_dens_m_km2
#min 0
#max 330306
Возможно, вы можете попробовать это еще раз в пустой папке. Я предполагаю, что в вашей папке остался дополнительный файл от предыдущей работы (например, Grip4_total_dens_m_km2.xml), содержащий эти числа. Возможно, этот файл был создан QGIS, возможно, с использованием быстрой оценки. Также сообщите, пожалуйста, версию terra, которую вы используете.
Спасибо за потрясающий пакет @RobertHijmans. Я сделал, как вы предлагаете, и получил те же результаты, что и вы. Я использовал терру 1.7.78. Извините, что отнимаю ваше время, но ценю ваш ответ. Я не уверен на 100%, что вызвало проблему, но думаю, что это может быть связано с сохранением растров в виде файлов .rda и последующей перезагрузкой в глобальную среду для использования. Я постоянно получал эту ошибку, которая, вероятно, вызывала проблему: Ошибка в .External(list(name = "CppMethod__invoke_notvoid", адрес = <pointer: (nil)>, : значение NULL, переданное как адрес символа. Спасибо за вашу помощь!
Вы можете сохранить SpatRaster с помощью writeRDS, но я не рекомендую это делать (вместо этого используйте .tif). Вы не можете использовать save (вероятно, это то, что вы имеете в виду под «.rda»); объекты будут повреждены, когда вы load их снова. Но это не связано с проблемой, о которой вы сообщили.
Это в некоторой степени связано с исходным вопросом: я столкнулся с проблемой при попытке упаковать различные наборы пространственных данных в пакет данных, который затем можно было бы загрузить и использовать при моделировании. Рекомендация (r-pkgs.org/data.html) — сохранять файлы данных как .rda (что дает сохранение). Путем проб и ошибок и ваших предложений я обнаружил, что сохранение растров в формате, отличном от .tif, проблематично. Это делает их упаковку намного сложнее, а доступ к ним с помощью пакета очень неуклюжим. Вряд ли стоит тратить усилия на упаковку. Спасибо за вашу помощь и отличные пакеты.
С помощью пакета вы можете хранить данные в любом формате, и написать функцию, которая возвращает их в желаемой форме, не составляет труда. Например: elev_data <- \() terra::rast(system.file("ex/elev.tif", package = "terra")).
некоторые сводные функции в R просматривают только случайную выборку данных для создания сводной статистики, если данные превышают определенный размер. посмотрите сюда