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

Исходный плакат рассчитывал примерно на 1 км на каждые 10 пикселей, но, вероятно, это слишком много. Я собираюсь создать сетку с размером ячейки 1 км * 1 км. Кроме того, исходный плакат не указывал источник сетки, поэтому я потрачу некоторое время на определение хорошей отправной точки. Я также предполагаю, что система координат проекции Сферический Меркатор является подходящим выбором для проекции. Это обычная проекция для Google Map или Open Street Maps.

1. Загрузить пакеты

Я собираюсь использовать следующие пакеты. Пакеты sp, rgdal и raster предоставляют множество полезных функций для пространственного анализа. leaflet и mapview - это пакеты для быстрой исследовательской визуализации пространственных данных.

# Load packages
library(sp)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)

2. Исследовательская визуализация местоположения станций.

Я создал интерактивную карту, чтобы осмотреть расположение четырех станций. Поскольку на исходном плакате указаны широта и долгота этих четырех станций, я могу создать SpatialPointsDataFrame с проекцией Широта / Долгота. Обратите внимание, что код EPSG для проекции широты / долготы - 4326. Чтобы узнать больше о коде EPSG, см. Это руководство (https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf).

# Create a data frame showing the **Latitude/Longitude**
station <- data.frame(lat = c(7.16, 8.6, 8.43, 8.15),
                      long = c(124.21, 123.35, 124.28, 125.08),
                      station = 1:4)

# Convert to SpatialPointsDataFrame
coordinates(station) <- ~long + lat

# Set the projection. They were latitude and longitude, so use WGS84 long-lat projection
proj4string(station) <- CRS("+init=epsg:4326")

# View the station location using the mapview function
mapview(station)

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

3. Определите происхождение

Изучив карту, я решил, что начало координат может быть около долготы 123 и широты 7. Это начало будет в нижнем левом углу сетки. Теперь мне нужно найти координату, представляющую ту же точку в сферической проекции Меркатора.

# Set the origin
ori <- SpatialPoints(cbind(123, 7), proj4string =  CRS("+init=epsg:4326")) 
# Convert the projection of ori
# Use EPSG: 3857 (Spherical Mercator)
ori_t <- spTransform(ori, CRSobj = CRS("+init=epsg:3857"))

Сначала я создал объект SpatialPoints на основе широты и долготы источника. После этого я использовал spTransform для преобразования проекта. Объект ori_t теперь является исходной точкой со сферической проекцией Меркатора. Обратите внимание, что код EPSG для сферического Меркатора - 3857.

Чтобы увидеть значение координат, мы можем использовать функцию coordinates следующим образом.

coordinates(ori_t)
     coords.x1 coords.x2
[1,]  13692297  781182.2

4. Определите размер сетки.

Теперь мне нужно определить размер сетки, которая может охватывать все четыре точки, и желаемую область для кригинга, которая зависит от размера ячейки и количества ячеек. Следующий код устанавливает экстент на основе информации. Я решил, что размер ячейки составляет 1 км * 1 км, но мне нужно поэкспериментировать с тем, какое число ячеек будет хорошим как для x-, так и для y-направления.

# The origin has been rounded to the nearest 100
x_ori <- round(coordinates(ori_t)[1, 1]/100) * 100
y_ori <- round(coordinates(ori_t)[1, 2]/100) * 100

# Define how many cells for x and y axis
x_cell <- 250
y_cell <- 200

# Define the resolution to be 1000 meters
cell_size <- 1000

# Create the extent
ext <- extent(x_ori, x_ori + (x_cell * cell_size), y_ori, y_ori + (y_cell * cell_size)) 

На основе созданного мной экстента я могу создать растровый слой с номерами, равными 0. Затем я могу снова использовать функцию mapview, чтобы проверить, хорошо ли совпадают растр и четыре станции.

# Initialize a raster layer
ras <- raster(ext)

# Set the resolution to be
res(ras) <- c(cell_size, cell_size)
ras[] <- 0

# Project the raster
projection(ras) <- CRS("+init=epsg:3857")

# Create interactive map
mapview(station) + mapview(ras)

Я повторил этот процесс несколько раз. В конце концов я решил, что количество ячеек составляет 250 и 200 для направлений x и y соответственно.

5. Создайте пространственную сетку.

Теперь я создал растровый слой нужного размера. Сначала я могу сохранить этот растр как GeoTiff для будущего использования.

# Save the raster layer
writeRaster(ras, filename = "ras.tif", format="GTiff") 

Наконец, чтобы использовать функции кригинга из пакета gstat, мне нужно преобразовать растр в SpatialPixels.

# Convert to spatial pixel
st_grid <- rasterToPoints(ras, spatial = TRUE)
gridded(st_grid) <- TRUE
st_grid <- as(st_grid, "SpatialPixels")

st_grid - это SpatialPixels, который можно использовать в кригинге.

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

@yzw и @Edzer помогают создать правильную прямоугольную сетку, но иногда возникает необходимость создать нерегулярную сетку поверх определенного многоугольника, обычно для кригинга. .

Это плохо документированная тема. Здесь можно найти один хороший ответ. Я расширяю его с помощью кода ниже:

Рассмотрим встроенный набор данных meuse. meuse.grid представляет собой сетку неправильной формы. Как сделать сетку вроде meuse.grid для нашей уникальной области изучения?

library(sp)
data(meuse.grid)
ggplot(data = meuse.grid) + geom_point(aes(x, y))

Представьте себе SpatialPolygon или SpatialPolygonsDataFrame неправильной формы, называемый spdf. Сначала вы строите над ним регулярную прямоугольную сетку, а затем подмножество точек в этой регулярной сетке с помощью многоугольника неправильной формы.

# First, make a rectangular grid over your `SpatialPolygonsDataFrame`
grd <- makegrid(spdf, n = 100)
colnames(grd) <- c("x", "y")

# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
  coords      = grd, 
  proj4string = CRS(proj4string(spdf))
)

# subset all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[spdf, ]


# Then, visualize your clipped grid which can be used for kriging
ggplot(as.data.frame(coordinates(grd_pts_in))) +
  geom_point(aes(x, y))

Если у вас есть область исследования в виде многоугольника, импортированного как SpatialPolygons, вы можете либо использовать пакетный растр для его растеризации, либо использовать sp::spsample для его выборки с использованием типа выборки regular.

Если у вас нет такого многоугольника, вы можете создавать точки, равномерно распределенные по прямоугольной области долготы / широты, используя expand.grid, используя seq для генерации последовательности длинных и широтных значений.