У меня есть климатические данные EURO-CORDEX, которые находятся на 11-градусной поворотной сетке полюсов. Я предварительно подготовил эти данные, преобразовав прогноз в WGS84. Данные поступают в виде точек, которые представляют центроиды квадратной сетки. Мне нужно создать квадратную сетку, которая обходит эти точки. Я получил общий метод для достижения этой цели, но конечные области ячеек сетки показывают ошибки до 50%.
Мой код указан ниже. Ранее я был предупрежден за предоставление кода в нотации Tidyverse, поэтому я стремился исключить это, где это возможно. Данные и код находятся на github здесь: https://github.com/avisserquinn/exampleData
Сначала я загружаю по долготе и широте моих центроидов из CSV и преобразовываю их во фрейм данных пространственных объектов с проекцией WGS84. Эти точки должны представлять сетку 11 x 11 градусов или 12 x 12 км.
> library(tidyverse)
> library(sf)
>
> data <- read_csv("stackExample.csv", col_types = cols())
> data <- st_as_sf(data, coords = c("lon", "lat")) # Spatial feature data frame
> data <- st_set_crs(data, 4326) # Set projection
> data
Simple feature collection with 2221 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: -9.996 ymin: 50.051 xmax: 1.965 ymax: 61.938
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
# A tibble: 2,221 x 2
grid geometry
<dbl> <POINT [°]>
1 1 (-9.996 51.768)
2 2 (-9.979 53.544)
3 3 (-9.96 52.013)
4 4 (-9.931 51.666)
5 5 (-9.924 52.258)
6 6 (-9.912 53.442)
7 7 (-9.906 54.034)
8 8 (-9.895 51.91)
9 9 (-9.875 53.687)
10 10 (-9.869 54.278)
# ... with 2,211 more rows
>
> ggplot(data) + geom_sf() + theme_bw()
Я создаю сеткуприменяя st_make_grid дважды (из пакета пространственных объектов sf). Первый раз нахожу центры между точками. Во второй раз я нахожу углы сетки, такие, что точки теперь являются центроидами.
> cellsize = .11
> dataGrid <- st_make_grid(data, cellsize = cellsize, what = "centers")
> dataGrid <- st_make_grid(dataGrid, cellsize = cellsize)
> dataGrid <- dataGrid %>% as_tibble %>% st_as_sf
> dataGrid
Simple feature collection with 11772 features and 0 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: -9.941 ymin: 50.106 xmax: 1.939 ymax: 62.096
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
First 10 features:
geometry
1 POLYGON ((-9.941 50.106, -9...
2 POLYGON ((-9.831 50.106, -9...
3 POLYGON ((-9.721 50.106, -9...
4 POLYGON ((-9.611 50.106, -9...
5 POLYGON ((-9.501 50.106, -9...
6 POLYGON ((-9.391 50.106, -9...
7 POLYGON ((-9.281 50.106, -9...
8 POLYGON ((-9.171 50.106, -9...
9 POLYGON ((-9.061 50.106, -8...
10 POLYGON ((-8.951 50.106, -8...
Затем я объединяю эти данные в сетке с центроидами, чтобы найти только подходящую сеткуЯчейки.
> dataGrid <- aggregate(data, dataGrid, FUN = mean)
> dataGrid <- as_tibble(dataGrid)
> dataGrid <- dataGrid[!is.na(dataGrid$grid),]
> dataGrid$area_sqm = st_area(dataGrid$geometry)
> dataGrid$area_sqkm = as.numeric(unlist(dataGrid$area_sqm * 10^-6))
> dataGrid$area_deficit = (12*12) - dataGrid$area_sqkm
> dataGrid
# A tibble: 2,175 x 5
grid geometry area_sqm area_sqkm area_deficit
<dbl> <POLYGON [°]> [m^2] <dbl> <dbl>
1 656 ((-5.651 50.106, -5.541 50.106, -5.541 50.216, -5.651 50.216, -5.651 50.106)) 96173304 96.2 47.8
2 678 ((-5.431 50.106, -5.321 50.106, -5.321 50.216, -5.431 50.216, -5.431 50.106)) 96173304 96.2 47.8
3 702 ((-5.211 50.106, -5.101 50.106, -5.101 50.216, -5.211 50.216, -5.211 50.106)) 96173304 96.2 47.8
4 730 ((-5.101 50.106, -4.991 50.106, -4.991 50.216, -5.101 50.216, -5.101 50.106)) 96173304 96.2 47.8
5 693 ((-5.321 50.216, -5.211 50.216, -5.211 50.326, -5.321 50.326, -5.321 50.216)) 95954257 96.0 48.0
6 720 ((-5.101 50.216, -4.991 50.216, -4.991 50.326, -5.101 50.326, -5.101 50.216)) 95954257 96.0 48.0
7 762 ((-4.881 50.216, -4.771 50.216, -4.771 50.326, -4.881 50.326, -4.881 50.216)) 95954257 96.0 48.0
8 1044 ((-3.891 50.216, -3.781 50.216, -3.781 50.326, -3.891 50.326, -3.891 50.216)) 95954257 96.0 48.0
9 712 ((-5.211 50.326, -5.101 50.326, -5.101 50.436, -5.211 50.436, -5.211 50.326)) 95734844 95.7 48.3
10 746. ((-4.991 50.326, -4.881 50.326, -4.881 50.436, -4.991 50.436, -4.991 50.326)) 95734844 95.7 48.3
# ... with 2,165 more rows
Когда я строю окончательный вывод, вы можете увидеть проблему. Размеры ячеек сетки значительно отклонены (дефицит);Я ожидаю некоторого отклонения от 12 х 12 км, так как проекция изогнута, но этот уровень дефицита крайне. Кроме того, в сетке есть пробелы;Я предполагаю, что это потому, что неправильный размер ячеек сетки означает, что не все точки захватываются?
> ggplot(dataGrid) +
+ geom_sf(aes(geometry = geometry, fill = area_deficit)) +
+ theme_minimal() +
+ scale_fill_viridis_c(direction = -1) +
+ scale_colour_viridis_c(direction = -1) +
+ coord_sf()
IЯ пробовал несколько вариантов, чтобы исправить это безуспешно. Будем весьма благодарны за любые советы по альтернативным решениям.