У меня есть очень простой вопрос, пожалуйста. Мне нужно обрезать несколько растров MODIS LAI, используя SpatVector
. Используя приведенный ниже код, обрезанный растр и SpatVector
не совпадают, поэтому я предполагаю, что обрезка выполнена неправильно. Есть идеи? Большое спасибо.
library(terra)
library(geodata)
#get spatVector
es <- geodata::gadm(country = "ESP", level=1, path = ".")
catal <- es[es$NAME_1 == "Cataluña", ]
#get MODIS data
modis <- terra::rast("see link below")
modis <- modis*0.1
# terra::plot(modis)
modis_crs <- terra::crs(modis) #get crs of MODIS data
catal_2 <- terra::project(catal, modis_crs) #give the spatVector the same crs as MODIS
# terra::crs(catal_2)==terra::crs(modis) #check if the two crs are the same
modis_cr <- terra::crop(modis, catal_2) #crop (here without mask)
terra::plot(modis_cr) #plot cropped MODIS data
terra::plot(catal_2, add=TRUE, lwd=2) #plot spatVector with the MODIS crs
Растр MODIS LAI для запуска этого примера доступен здесь <- ссылка удалена
Вот что я получаю: два файла не совпадают (игнорируйте цвет растра):
Я пытался изменить crs данных MODIS, а не crs spatVector
, я пробовал использовать разные варианты extent
, align
или snap
. Я попробовал немного поиграться с библиотекой sf
и прочитал немало подобных вопросов.
Похоже, определение вашего растра в crs неверно. Переопределение этого значения с помощью приведенной ниже строки proj4string приводит к ожидаемому позиционированию:
library(terra)
#> terra 1.7.71
library(geodata)
# get data
es <- geodata::gadm(country = "ESP", level=1, path = ".")
catal <- es[es$NAME_1 == "Cataluña", ]
# read modis
modis <- terra::rast("cropmodis.tif")
# replace modis crs
terra::crs(modis) <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"
# reproject vector
modis_crs <- terra::crs(modis)
catal_2 <- terra::project(catal, modis_crs)
# crop
modis_cr <- terra::crop(modis, catal_2)
# inspect
terra::plot(modis_cr)
terra::plot(catal_2, add = TRUE)
Примечание:
Определение crs взято из случайного набора данных MODIS, загруженного с EarthExplorer:
terra::rast("MOD11A2.A2024089.h18v04.061.2024099144843.hdf") |> terra::crs(proj = TRUE)
#> [1] "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"
@Chris: Не обязательно, но количество символов здесь ограничено, чтобы опубликовать определение WKT2 ;-)
Это делает ?terra::crs
интересным чтение в контексте представленной проблемы. Допускает ли атрикуляция terra::crs использование wkt, ухудшенного в настоящее время. Во что-то поиграть.
Фантастика, это решило проблему. Большое спасибо. Эти данные MODIS LAI были загружены из Google Earth Engine, и поэтому я предполагаю, что GEE как бы потеряла правильную информацию CRS при объединении плиток для покрытия Испании. Еще раз большое спасибо.
@Chris: По крайней мере, WKT2, похоже, является представлением crs по умолчанию при использовании terra::crs()
, несмотря на подробное примечание в ?terra::crs
.
Ну, если по умолчанию, я бы описал это как reluctant default
: действуйте осторожно. Я думаю, что ваше и @Fran предположение об обработке GEE для объединенных плиток, скорее всего, является виновником. Что касается proj и WKT2, они кажутся независимыми, т.е. +proj не является производным от элементов W, поскольку GEE WKT2 имеет эллипсоид 6378137,298.257223563, который не втягивается в +R, как ваш Earth Explorer, а только + данные. Что ж, мы все будем знать, что касается будущих вопросов GEE.
@Chris: Я бы сказал, что это просто неправильное определение crs, и я не ожидал бы, что речь идет о типе (поскольку мы ничего не знаем о выполняемой предварительной обработке растра). Но, похоже, это какая-то операция слияния, чтобы покрыть территорию Испании отдельными фрагментами и преобразованием формата, поскольку данные MODIS не поступают в формате GeoTIFF, по крайней мере, на самом деле.
terra::crs(modis, proj = TRUE)
возвращается"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
кстати.