Блог

Я пытаюсь перепроектировать глобальный растр искусственного излучения света с долготы / широты на равную площадь Берманна (EPSG: 6933) с разрешением 0,0417 градусов. Из-за скачков данных в городских районах, когда пиксели интерполируются во время перепроектирования, происходит потеря данных на уровне около 15% по всему слою.

Я попытался преобразовать растр в фрейм данных пространственных точек, перепроектировать фрейм данных пространственных точек, а затем растрировать с использованием растра, созданного с использованием функции «projectraster» в качестве растра шаблона (я думаю, что размеры, протяженность и разрешение растра шаблона могут быть проблемой?) Однако при этом получается растр с горизонтальными линиями через слой.

Вот пример кода с Испанией в качестве примера. Я могу отправить по электронной почте файл tif для Испании (246 КБ):

 library("sf")
library("raster")

behrmann <- CRS(' proj=cea  lon_0=0  lat_ts=30  x_0=0  y_0=0  datum=WGS84  ellps=WGS84  units=m  no_defs')

r <- raster("~/Documents/R spatial data[enter image description here][1]/Spain.tif")
cellStats(r, sum) # check summed light emissions
r_temp <- projectRaster(r, crs = behrmann) # creates template for rasterisation (data is lost due to interpolation of data spikes)
spdf <- rasterToPoints(r, spatial = TRUE)
spdf2 <- spTransform(spdf, CRS = behrmann)
r2 <- rasterize(spdf2, r_temp, field = "Spain", fun = "sum")
cellStats(r2, sum) # check no data has been lost
plot(log10(r2)) # see attached image[enter image description here][1]
  

Как можно перепроектировать на равную область без потери данных и избегая горизонтальных линий? Я также попытался преобразовать в фрейм пространственных полигонов вместо пространственных точек, и это не приводит к созданию линий, а вместо этого теряет данные, аналогичные функции ‘projectRaster’. Это должно быть распространенной проблемой, но я не могу найти никакой помощи в Интернете.

Заранее большое спасибо.

Пример горизонтальных линий после перепроектирования stack.imgur.com/IV0fZ.png

При преобразовании растра вычисляются новые значения ячеек. Обычно это делается путем усреднения, что приводит к менее экстремальным значениям.

 library(raster)
r <- raster(res=5)
set.seed(1)
values(r) <- runif(ncell(r))
r <- focal(r, w=matrix(1, 3, 3))
r <- focal(r, w=matrix(1, 3, 3))
r <- round(focal(r, w=matrix(1, 3, 3)))
r
#class      : RasterLayer 
#dimensions : 36, 72, 2592  (nrow, ncol, ncell)
#resolution : 5, 5  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        :  proj=longlat  datum=WGS84  no_defs 
#source     : memory
#names      : layer 
#values     : 229, 473  (min, max)


behrmann <- " proj=cea  lon_0=0  lat_ts=30  x_0=0  y_0=0  datum=WGS84  ellps=WGS84  units=m"
x <- projectRaster(r, crs=behrmann)
x
#class      : RasterLayer 
#dimensions : 27, 78, 2106  (nrow, ncol, ncell)
#resolution : 482000, 638000  (x, y)
#extent     : -18813530, 18782470, -8607770, 8618230  (xmin, xmax, ymin, ymax)
#crs        :  proj=cea  lat_ts=30  lon_0=0  x_0=0  y_0=0  datum=WGS84  units=m  no_defs 
#source     : memory
#names      : layer 
#values     : 233.0725, 471.7214  (min, max)
  

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

 z <- projectRaster(r, crs=behrmann, method="ngb")
z
#class      : RasterLayer 
# ...
#values     : 229, 473  (min, max)
  

или действительно используйте полигоны (если у вас не слишком много ячеек)

 p <- as(r, "SpatialPolygonsDataFrame")
y <- spTransform(p, behrmann)

y
#class       : SpatialPolygonsDataFrame 
#features    : 2160 
#extent      : -17367530, 17367530, -7089914, 7089914  (xmin, xmax, ymin, ymax)
#crs         :  proj=cea  lat_ts=30  lon_0=0  x_0=0  y_0=0  datum=WGS84  units=m  no_defs 
#variables   : 1
#names       : layer 
#min values  :   229 
#max values  :   473