Перепроектировать растр (от WGS84 до BNG) с большим разрешением

Я делаю достаточно простую операцию, репроецируя растр из WGS84 в British National Grid, но мне интересно узнать о некоторых результатах после репроецирования. Сумма получаемых растров весьма различна; это связано с разрешением и с тем, как работает билинейная интерполяция в 'projectRaster'?

Я начал с CSV, который имеет глобальные данные, охватывающие диапазон от -180 до 180 градусов по широте и долготе (только целочисленные значения), с некоторыми значениями z. это подмножество, преобразованное в растр проекции WGS84 и преобразованное в BNG (подмножество ниже):

x <- rep(c(-10:3), times = 10)
n <- 14
y <- rep(48:57, each=n)
z <- rnorm(n=140, mean=20, sd=5)
ind <- which(z %in% sample(z, 45))
z[ind]<-NA
df <- data.frame("x"=x,"y"=y,"value"=z)

bng <- '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs'
wgs84 <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'

coordinates(df) = ~x + y
ras = rasterFromXYZ(df, crs=wgs84)
cellStats(ras,sum)
plot(ras)

ras_bng <- projectRaster(ras,crs=bng)
plot(ras_bng)
cellStats(ras_bng,sum)

сумма идет с 1919 по 2625 (в моем случае), довольно много.

Это только из-за репроекции, создающей так много дополнительных ячеек вокруг "краев"? если бы мне пришлось перепроецировать в растр с гораздо меньшим разрешением (5 км), это значительно уменьшило бы разные суммы?

спасибо, с

> ras
class       : RasterLayer 
dimensions  : 10, 14, 140  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -10.5, 3.5, 47.5, 57.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 3.243918, 32.21532  (min, max)

> ras_bng
class       : RasterLayer 
dimensions  : 12, 19, 228  (nrow, ncol, ncell)
resolution  : 67900, 111000  (x, y)
extent      : -375677.8, 914422.2, -343553.2, 988446.8  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 
data source : in memory
names       : layer 
values      : 8.174968, 28.31331  (min, max)

1 ответ

Решение

Это потому что NA значения в "RAS".

Разница мала, когда вы вынимаете z[ind]<-NA, Это потому что projectRaster неявно использует "na.rm=TRUE"; возможно, нужен аргумент, чтобы это изменить. Это не так просто, как в других случаях, так как, например, только одна из ячеек, интерполированных из NA, в этом случае это, вероятно, должно быть вычислено.

На практике это редко можно найти столько значений NA, как правило, они ограничены только краями (суши).

Кстати WGS84 это данные, это не проекция. Используемая вами проекция будет представлять собой "долготу / широту" (также известную под другими именами), за исключением того, что это также не проекция, поскольку вся точка проекции состоит в переходе от таких угловых координат к плоским координатам. Таким образом, вы преобразуете одну систему координат (long/lat, WGS84) в другую (BNG).

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