Проблема при обрезке растровых данных с помощью spatVector

У меня есть очень простой вопрос, пожалуйста. Мне нужно обрезать несколько растров 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 и прочитал немало подобных вопросов.

Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
2
0
99
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Похоже, определение вашего растра в 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: Я бы сказал, что это просто неправильное определение 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" кстати.

dimfalk 20.04.2024 01:57

@Chris: Не обязательно, но количество символов здесь ограничено, чтобы опубликовать определение WKT2 ;-)

dimfalk 20.04.2024 02:03

Это делает ?terra::crs интересным чтение в контексте представленной проблемы. Допускает ли атрикуляция terra::crs использование wkt, ухудшенного в настоящее время. Во что-то поиграть.

Chris 20.04.2024 04:22

Фантастика, это решило проблему. Большое спасибо. Эти данные MODIS LAI были загружены из Google Earth Engine, и поэтому я предполагаю, что GEE как бы потеряла правильную информацию CRS при объединении плиток для покрытия Испании. Еще раз большое спасибо.

Fran 21.04.2024 15:16

@Chris: По крайней мере, WKT2, похоже, является представлением crs по умолчанию при использовании terra::crs(), несмотря на подробное примечание в ?terra::crs.

dimfalk 21.04.2024 22:22

Ну, если по умолчанию, я бы описал это как reluctant default: действуйте осторожно. Я думаю, что ваше и @Fran предположение об обработке GEE для объединенных плиток, скорее всего, является виновником. Что касается proj и WKT2, они кажутся независимыми, т.е. +proj не является производным от элементов W, поскольку GEE WKT2 имеет эллипсоид 6378137,298.257223563, который не втягивается в +R, как ваш Earth Explorer, а только + данные. Что ж, мы все будем знать, что касается будущих вопросов GEE.

Chris 22.04.2024 02:37

Другие вопросы по теме