Подсчет полигонов внутри растровых ячеек в R

Содержание

Вопрос или проблема
Ответ или решение
Подготовка данных
Подсчет полигонов в ячейках растровой сетки
Обработка ошибок
Отображение результата
Альтернативный способ подсчета с использованием cellFromXY

Вопрос или проблема

Я пытаюсь сделать что-то, что, как мне казалось, легко, а именно подсчитать количество многоугольников в каждой растрной ячейке. Я пробую оба решения здесь: R: Как я могу подсчитать, сколько точек в каждой ячейке моей сетки?, но ни одно из них не работает, и я буквально схожу с ума от всех используемых различных пакетов. Я хотел бы решение на terra, но не знаю, как этого добиться, и в данный момент открыт к любым пакетам. Вот некоторые примерные данные и мой код на данный момент:

# Загружаем пакеты
pacman::p_load(terra, raster, rnaturalearth, rts, sf)

# Загружаем примерные многоугольники
cntrs <- ne_countries(country = c("united kingdom", "argentina", "singapore"), type = "countries")

# Проецируем
cntrs <- project(vect(cntrs), "+proj=moll +lon_0=0 +x_0=0 +y_0=0")

# Создаем пустую растрную сетку
r <- terra::rast(xmin = -180, xmax = 180, ymin = -90, ymax = 90,
          res = 10,
          crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

r[] <- 0

### Первое решение из ссылки: растризация
terra::rasterize(cntrs, r, fun = "count")

# Выдает предупреждающее сообщение и пустой растр
Предупреждающее сообщение:
  Не удалось вычислить min/max, не найдено действительных пикселей при выборке. (Ошибка GDAL 1)

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

### Второе решение: cellFromXY
# Изменяем формат многоугольников
cntrs <- as(st_geometry(st_as_sf(cntrs)), "Spatial")
table(cellFromXY(r, cntrs))

# Выдает сообщение об ошибке:
Ошибка в (функция (классы, fdef, mtable)  : 
            невозможно найти унаследованный метод для функции ‘cellFromXY’ для сигнатуры ‘"SpatRaster", "SpatialPolygons"’

Хотя моя версия rts rts_1.1-14, такая же, как в мануале (https://cran.r-project.org/web/packages/rts/rts.pdf), где они пишут, что “объект” может быть SpatRaster…

# Пытался преобразовать SpatRaster в RasterLayer, как в примерах по ссылке
r <- raster::raster(r)
table(cellFromXY(r, cntrs))

# Выдает сообщение об ошибке:
Ошибка: Не совместимо с запрашиваемым типом: [тип=S4; цель=double].

В примере они используют SpatialPoints, но у меня многоугольники, так что я пытаюсь с SpatialPolygons, что тоже, похоже, не работает. Буду признателен за любую помощь.

Ответ или решение

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

Подготовка данных

Загрузка необходимых библиотек: Вы правильно подключили пакеты, однако для наших нужд нам понадобится только terra и sf, поскольку terra может работать непосредственно с векторными данными, представленными в формате sf.
```
library(terra)
library(sf)
```
Загрузка и преобразование полигонов: Используем функцию ne_countries для получения полигонов стран. При этом, если данные изначально в формате sf, можно обойтись без их преобразования.
```
cntrs <- ne_countries(country = c("united kingdom", "argentina", "singapore"), type = "countries", returnclass = "sf")
```
Пространственная проекция: Убедитесь, что ваши полигоны имеют необходимые координаты.
```
cntrs <- st_transform(cntrs, crs = "EPSG:3395")  # Используем проекцию Меркатора
```

Создание растровой сетки:

r <- rast(xmin = -180, xmax = 180, ymin = -90, ymax = 90, res = 10, crs = "EPSG:3395")

Подсчет полигонов в ячейках растровой сетки

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

# Подсчёт количества полигонов в ячейках растровой сетки
result <- terra::rasterize(cntrs, r, fun = "count", background = 0)

Обратите внимание на параметр background. Если полигоны не перекрывают ячейки, мы бы хотели установить значение по умолчанию в 0, чтобы отражать отсутствие полигонов в конкретной ячейке.

Обработка ошибок

Если вы всё ещё сталкиваетесь с проблемой "Failed to compute min/max, no valid pixels found in sampling", это может свидетельствовать о том, что ни один из полигонов не пересекает растровую ячейку. Также возможно, что CRS растровой сетки и полигонов не совпадают, поэтому всегда важно проверять это.

Отображение результата

Для визуализации результата, вы можете использовать функцию plot:

plot(result, main = "Количество полигонов в ячейках растровой сетки")

Альтернативный способ подсчета с использованием `cellFromXY`

Если вы хотите использовать метод cellFromXY, вам необходимо преобразовать полигоны в их центры и затем использовать их для подсчета:

# Получение центров полигонов
centroids <- st_centroid(cntrs)

# Получение ячеек растровой сетки для центров
cell_indices <- cellFromXY(r, st_coordinates(centroids))

# Подсчет количества полигонов в каждой ячейке
counts <- table(cell_indices)

# Отображение результатов
print(counts)

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