разные результаты при проецировании растра с агрегатом или без него - PullRequest
Новая
       26

разные результаты при проецировании растра с агрегатом или без него

0 голосов
/ 24 октября 2018

Моя цель - перепроектировать растр с именем pftnc, чтобы иметь тот же CRS и размер ячейки, что и у другого растра sp, чтобы они идеально перекрывались и могли быть сложены.Доступные здесь файлы pftnc и sp

Я думал, что это проблема с программным обеспечением, см. Другой вопрос https://gis.stackexchange.com/questions/299935/extract-values-from-raster-r-and-arcgis-different-results,, но у меня возникают несоответствия только в пределахR. Итак, я хотел бы сначала найти решение R.

Я попытался перепроектировать напрямую 

pft1 <- projectRaster(pftnc, sp)

А также попытался сначала агрегировать, а затем перепроецировать 

pftagg <- aggregate(pftnc, fact=4)
pftproj1 <- projectRaster(from=pftagg, to=sp)

Я также пытался просто спроецировать, а затем обрезать и агрегировать, или наоборот, но я никогда не смог получить результат, который получаю, когда использую to = sp, поэтому я отказался от этого маршрута.Например, 

cea=crs(sp)
pftproj2 <- projectRaster(pftnc, crs=cea)

Проблема заключается в том, что pftproj1 и pft1 имеют одинаковые строки, столбцы, экстент и разрешение, но значения отличаются. 

> cellStats(pftproj1,mean)
[1] 0.1955
> cellStats(pft1,mean)
[1] 0.2028

Существует ли правильный порядокоперации для достижения этой цели?Должен ли я позволить projectRaster иметь дело с агрегацией значений?Я представляю

1 Ответ

0 голосов
/ 31 октября 2018

Я не знаю, есть ли «лучший» метод.Я думаю, что основная причина того различия, которое вы получаете, состоит в том, что если вы сначала агрегируете, вы получаете больше ячеек (если агрегируете с na.rm=TRUE).Разница намного меньше, если вы используете mask для сравнения строго одинаковых ячеек. 

library(raster)
sp <- raster("sp.tif")
pf <- raster("pftnc.tif")

pf1 <- projectRaster(pf, sp)
pfa <- aggregate(pf, fact=4)
pf2 <- projectRaster(from=pfa, to=sp)
pf2m <- mask(pf2, pf1)

# an alternative that I would not recommend
pfr <- projectRaster(pf, crs=crs(sp))
pf3 <- resample(pfr, sp)
pf3m <- mask(pf3, pf1)

cellStats(pf, mean)
#[1] 0.2202417
cellStats(pf1, mean)
#[1] 0.2027533
cellStats(pf2, mean)
#[1] 0.1954503
cellStats(pf2m, mean)
#[1] 0.2024068
cellStats(pf3, mean)
#[1] 0.2016164
cellStats(pf3m, mean)
#[1] 0.2021608

В вопросе лучше использовать пример данных, сгенерированных кодом.В вашем случае вы могли бы использовать 

sp <- raster(nrow=142, ncol=360, ext=extent(-17367529, 17367529, -6356742, 7348382), crs="+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m")

pf <- raster(ncol=1440, nrow=720, xmn=-180, xmx=180, ymn=-90, ymx=90, crs='+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
values(pf) = rep(c(1:4,NA), ncell(pf)/5)

Но ... различия исчезли, за исключением «плохого» метода.Не уверен, что это связано с распределением ваших значений или с ролью значений NA, но я думаю, вы могли бы поэкспериментировать, чтобы выяснить это.

cellStats(pf, mean)
#[1] 2.5
cellStats(pf1, mean)
#[1] 2.5
cellStats(pf2, mean)
#[1] 2.5
cellStats(pf2m, mean)
#[1] 2.5
cellStats(pf3, mean)
#[1] 2.547783
cellStats(pf3m, mean)
[1] 2.547783
...