Создать Grid в R для кригинга в GSTAT
lat long
7.16 124.21
8.6 123.35
8.43 124.28
8.15 125.08
Рассмотрим эти координаты, эти координаты соответствуют метеостанциям, которые измеряют данные об осадках.
Вступление к пакету gstat в R использует набор данных meuse. В какой-то момент этого урока: https://rpubs.com/nabilabd/118172, ребята используют "meuse.grid" в этой строке кода:
data("meuse.grid")
У меня нет такого файла, и я не знаю, как его создать, могу ли я создать его, используя эти координаты? Или, по крайней мере, укажите мне материал, который обсуждает, как создать пользовательскую сетку для пользовательской области (то есть, не используя административные границы из GADM).
Вероятно, это неправильно, даже не знаю, имеет ли этот вопрос смысл для здравомыслящих людей. Тем не менее, хотелось бы услышать какое-то направление или хотя бы советы. Большое спасибо!
Итого нуб на R и статистика.
РЕДАКТИРОВАТЬ: Посмотрите пример сетки, как выглядит учебник, который я опубликовал, вот что я хочу сделать.
РЕДАКТИРОВАТЬ 2: этот метод будет жизнеспособным? https://rstudio-pubs-static.s3.amazonaws.com/46259_d328295794034414944deea60552a942.html
3 ответа
Я собираюсь поделиться своим подходом к созданию сетки для кригинга. Вероятно, есть более эффективные или элегантные способы решения той же задачи, но я надеюсь, что это станет началом некоторых дискуссий.
Оригинальный плакат думал о 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=x, y=y))
Вообразите неправильную форму SpatialPolygon или же SpatialPolygonsDataFrame, называется spdf, Сначала вы строите правильную прямоугольную сетку поверх нее, а затем подставляете точки в этой регулярной сетке многоугольником неправильной формы.
1. Сделайте прямоугольную сетку над вашим SpatialPolygonsDataFrame
grd <- makegrid(spdf, n = 100)
colnames(grd) <- c('x','y')
2. Преобразовать сетку в SpatialPoints и подстановить эти точки многоугольником.
grd_pts <- SpatialPoints(coords = grd,
proj4string=proj4string(spdf))
# find all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[spdf, ]
3. Визуализируйте свою обрезанную сетку, которую можно использовать для кригинга.
# transform grd_pts_in back into a data frame
gdf <- as.data.frame(coordinates(grd_pts_in))
ggplot(gdf) +
geom_point(aes(x=x,y=y))
Если у вас есть учебная область в виде многоугольника, импортируется как SpatialPolygonsВы можете использовать растровый пакет для его растеризации или использовать sp::spsample пробовать его, используя тип выборки regular,
Если у вас нет такого многоугольника, вы можете создавать точки, регулярно растягивающиеся по прямоугольной области длинных / лат, используя expand.grid, с помощью seq генерировать последовательность значений long и lat.