Для функции я хочу перепроектировать растровый ввод - сам вывод после использования кадрирования и маски - с использованием установленного пользователем CRS. Я думал, что проецирование с существующими crs ничего не сделает и просто вернет входной растр. К моему удивлению, это был не тот случай. Ниже воспроизводимого примера
Создание фиктивного растра:
library(raster)
# Download country polygin, in this case Malawi
mwi <- raster::getData("GADM", country = "MWI", level = 0)
# Create dummy raster
grid <- raster::raster() # 1 degree raster
grid <- raster::disaggregate(grid, fact = 12) # 5 arcmin
grid <- raster::crop(grid, mwi)
values(grid) <- rep(1, ncell(grid)) # Set a value
# The input raster with dimensions 94,39,3666
grid <- raster::mask(grid, mwi)
plot(grid)
grid
# Reproject the raster using its own crs. I use ngb as it is a categorical variable.
# This raster has dimensions 102, 47, 4794 so it seems a lot of white space (NA) is added.
own_crs <- crs(grid)
grid_reproj <- raster::projectRaster(grid, crs = own_crs, method = "ngb")
plot(grid_reproj)
grid_reproj
# To remove the white space I use trim
# This results in a waster with dimensions 93, 39, 3627
grid_trim <- raster::trim(grid_reproj)
plot(grid_trim)
grid_trim
# I also decided to compare the maps visually with mapview
library(mapview)
# There seems to be a trim function in mapview which I set to FALSE
# Also use the browser for easy viewing
options(viewer = NULL)
mapviewOptions(trim = FALSE)
mapviewGetOption("trim")
mapview(grid, col.regions = "green", legend = F) +
mapview(grid_reproj, col.regions = "red", legend = F) +
mapview(grid_trim, col.regions = "blue", legend = F)
При сравнении карт я наблюдаю две вещи:
(1) grid и grid_trim почти идентичны, за исключением одного ячейка сетки сверху. Возможно, это связано с округлением,
(2) grid_reproj имеет гораздо большие размеры и разную степень. Также кажется, что карта немного смещена по сравнению с другими картами. Это корректируется с помощью обрезки, поэтому я предполагаю, что это на самом деле ячейки NA, и разница может быть связана с тем, как mapview отображает карты.
Следовательно, мой главный вопрос заключается в том, что происходит, когда rasterProject проектирует с тем же степени? И почему результат отличается даже после обрезки?