Растровое изображение кажется смещенным при использовании листочка для R

Я хочу построить некоторые пространственные данные с помощью пакета листовок в R, однако сгенерированное растровое изображение кажется смещенным по сравнению с опорной сеткой. Я подозреваю, что проблема с проекцией карты, но я не разбираюсь в этой теме, поэтому любая помощь будет принята с благодарностью.

Вот минимальный код для построения карты:

library(leaflet)
library(sp)
library(raster)

set.seed(111)

# create dummy data -rectangular grid with random values
m = 10
n = 10
x = seq(45,48,length.out = m)
y = seq(15,18,length.out = n)
X = matrix(rep(x, each = n), nrow = n)
Y = matrix(rep(y, m), nrow = n)

# collector dataframe
points = data.frame(value = rnorm(n*m), lng = c(Y), lat = c(X))

## create raster grid
s = SpatialPixelsDataFrame(points[,c('lng', 'lat')], data = points)
# set WGS84 projection
crs(s) = sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

r = raster(s)

# add coloring
pal = colorNumeric(c("#0C2C84", "#f7f7f7", "#F98009"), points$value,
                    na.color = "transparent")

## plot map
leaflet() %>% addProviderTiles("CartoDB.Positron") %>%
   addRasterImage(r, colors = pal, opacity = 0.6)

В результате получается следующая карта, которая на первый взгляд подходит: введите описание изображения здесь

Если на эту карту добавлена ​​сетка:

## grid
dx = diff(x)[1]/2
dy = diff(y)[1]/2

rect_lng = sapply(points$lng, function(t) c(t-dx, t+dx))
rect_lat = sapply(points$lat, function(t) c(t-dy, t+dy))

leaflet() %>% addProviderTiles("CartoDB.Positron") %>% 
  addRectangles(
    lng1=rect_lng[1,], lat1=rect_lat[1,],
    lng2=rect_lng[2,], lat2=rect_lat[2,],
    fillColor = "transparent",
    weight = 1
  ) %>%
  addRasterImage(r, colors = pal, opacity = 0.6)

Карта выглядит так:

введите здесь описание изображения

Здесь мы видим, что сетки не совпадают. введите здесь описание изображения

В чем причина этого несоответствия? Как это можно было устранить? Зря пробовал разные проекции. Единственное, что сработало, - это использовать addRectangle вместо addRasterImage, однако это требует гораздо больше вычислений и замедляет процесс, поэтому я хочу избежать. Обратите внимание, что в приведенном выше примере addRectangle используется только для ссылки, в окончательном коде я не хочу его использовать.

Для карт с большим количеством ячеек (сеток) несоответствие довольно велико, может превышать размер одной ячейки.

Заранее благодарю за любую помощь.

ИЗМЕНИТЬ

Проблема может быть связана с проблемой проекции между проекциями эллипсоида и сферы, см. Последний вопрос здесь :

для преобразования между WGS84 и меркатором на сфере будут существенные сдвиги в координатах Y-меркатора. Это связано с тем, что внутренне cs2cs должен настроить координаты широты и долготы с точки на сфере до точки привязки WGS84, которая имеет эллипсоид совершенно другой формы.

Однако я не смог решить проблему с помощью рекомендованного «трюка»: +nadgrids=@null.


r leaflet sp r-raster map-projections
person rozsasarpi    schedule 16.11.2015    source источник

Ответы (1)


arrow_upward
4
arrow_downward

Автор буклета R здесь. Мне кажется, что средство визуализации растровых слоев, которое я написал, начинает дрейфовать, когда исходный растр имеет очень мало пикселей по сравнению с количеством пикселей, отображаемых на экране. Вы можете убедиться в этом, внеся следующие изменения в растр:

r1 <- r
nrow(r1) <- 600
ncol(r1) <- 600
r <- resample(r, r1, method = "ngb")

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

person Joe Cheng    schedule 17.11.2015
comment
Я думаю, проблема в том, что выполняемая мной операция проецирования (leaflet::projectRasterForLeaflet) создает выходной растр с теми же размерами пикселей / ячеек, что и входной растр (в данном случае 10x10), тогда как обычно вы хотите проецировать с использованием выходных размеров, которыми вы на самом деле являетесь. рендеринг в. Но в этом случае на самом деле нет фиксированных выходных размеров, поскольку пользователь может увеличивать и уменьшать масштаб. Я делаю небольшой проецируемый растр, а затем увеличиваю его линейно, чтобы соответствовать зуму, и, вероятно, это то, что вводит дрейф. - person Joe Cheng; 17.11.2015
comment
Спасибо за обходной путь и возможное объяснение! Растр с повторной дискретизацией r ведет себя странно в colorNumeric; Я изменил домен на values(r), но получаю предупреждение о том, что некоторые значения находятся за пределами цветовой шкалы и действительно некоторые пиксели прозрачны (дыры на растровом графике). В приведенном выше примере domain = c(1.01*min(values(r)), 1.01*max(values(r))) решил проблему, но более элегантное решение было бы лучше. Если у вас есть предложения по этому поводу, это было бы здорово. Спасибо! - person rozsasarpi; 17.11.2015
comment
Хм. domain = c(minValue(r), maxValue(r)) не работает (где r - растр с уже передискретизацией)? Зачем нужен 1.01? - person Joe Cheng; 17.11.2015
comment
Нет, domain = c(minValue(r), maxValue(r)) не работает, поэтому я применил множитель 1.01 (пробовал с передискретизированным r). - person rozsasarpi; 17.11.2015