Boxplot растра и шейп-файла
Этот пример кода работает
library(raster)
library(sp)
library(rgdal)
library(ggplot2)
# Create some sample raster data
raster_file <- raster(ncol=36, nrow=18)
raster_file[] <- 1:ncell(raster_file)
plot(raster_file)
#Create some sample polygons
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
shape_file <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1),
Polygons(list(Polygon(cds2)), 2)))
plot(shape_file)
# Extract raster values within the shapefile extracted_values <- extract(raster_file, shape_file)
# Assuming the shapefile has multiple polygons and you want to create a boxplot for each
data_list <- lapply(1:length(extracted_values), function(i) {
data.frame(value = extracted_values[[i]], polygon = i)
})
data <- do.call(rbind, data_list)
# Create the boxplot
bp<-ggplot(data, aes(x = factor(polygon), y = value)) +
geom_boxplot() +
labs(x = "Polygon", y = "Raster Values") +
theme_minimal()
bp
В моем собственном наборе данных я столкнулся с проблемами чтения полигонов.
Сообщение об ошибке появляется в функции boxplot
Здесь может быть файл формы
> # load shape file including all layers and print layers
> shape_file<-shapefile("C:/Users/.....BiogeoRegion.shp")
Warning message:
[vect] Z coordinates ignored
> names(shape_file)
[1] "RegionNumm" "RegionName" "Unterregio" "Unterreg_1" "ObjNummer"
"Version" "Shape_Leng" "Shape_Area" "DERegionNa" "FRRegionNa"
[11] "ITRegionNa" "DEBioBedeu" "FRBioBedeu" "ITBioBedeu"
> str(shape_file)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
..@ data :'data.frame': 12 obs. of 14 variables:
.. ..$ RegionNumm: int [1:12] 1 2 2 2 2 3 3 4 5 5 ...
.. ..$ RegionName: chr [1:12] "R1" "R2" "R2" "R2" ...
.. ..$ Unterregio: int [1:12] 11 21 22 23 24 31 32 41 51 52 ...
.. ..$ Unterreg_1: chr [1:12] "U11" "U21" "U22" "U23" ...
.. ..$ ObjNummer : chr [1:12] "1" "2" "3" "4" ...
.. ..$ Version : chr [1:12] "2020/05/08" "2020/05/08" "2020/05/08"
"2020/05/08" ...
> # select a specific raster file from a list
> raster_file<-allrasters_pres[[1]]
> names(raster_file)
[1] "Andrena.barbilabris_glo_ensemble"
>
> # Extract raster values within the shapefile extracted_values <-
> extract(raster_file, shape_file)
>
>
>
> # Assuming the shapefile has multiple polygons and you want to create
> a
boxplot for each
> data_list <- lapply(1:length(extracted_values), function(i) {
+ data.frame(value = extracted_values[[i]], polygon = i)
+ })
> data <- do.call(rbind, data_list)
> names(data)
[1] "value" "polygon"
>
> # Create the boxplot
> bp<-ggplot(data, aes(x = factor(polygon), y = value)) +
+ geom_boxplot() +
+ labs(x = "Polygon", y = "Raster Values") +
+ theme_minimal()
> bp
Error in UseMethod("depth") :
no applicable method for 'depth' applied to an object of class "NULL"
In addition: Warning message:
Removed 452451 rows containing non-finite outside the scale range (`stat_boxplot()`).
>
Ответы 2
Если это не обязательно должно быть основано на пакетах raster и sp, вы можете вместо этого рассмотреть terra или sf + stars.
С помощью копипаста из документации terra это может выглядеть примерно так:
library(terra)
#> terra 1.7.78
library(ggplot2)
# Load shp and raster from example files of terra
(v <- vect(system.file("ex/lux.shp", package = "terra")))
#> class : SpatVector
#> geometry : polygons
#> dimensions : 12, 6 (geometries, attributes)
#> extent : 5.74414, 6.528252, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
#> source : lux.shp
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> names : ID_1 NAME_1 ID_2 NAME_2 AREA POP
#> type : <num> <chr> <num> <chr> <num> <int>
#> values : 1 Diekirch 1 Clervaux 312 18081
#> 1 Diekirch 2 Diekirch 218 32543
#> 1 Diekirch 3 Redange 259 18664
(r <- rast(system.file("ex/elev.tif", package = "terra")))
#> class : SpatRaster
#> dimensions : 90, 95, 1 (nrow, ncol, nlyr)
#> resolution : 0.008333333, 0.008333333 (x, y)
#> extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source : elev.tif
#> name : elevation
#> min value : 141
#> max value : 547
plot(r); plot(v, add = TRUE); text(v, labels = "NAME_2", col=c("white"))
# extract, include names from vector attributes
e <- extract(r,v)
e$NAME_2 <- v$NAME_2[e$ID]
str(e)
#> 'data.frame': 4606 obs. of 3 variables:
#> $ ID : num 1 1 1 1 1 1 1 1 1 1 ...
#> $ elevation: int 547 485 497 515 515 515 518 531 540 NA ...
#> $ NAME_2 : chr "Clervaux" "Clervaux" "Clervaux" "Clervaux" ...
# plot
ggplot(e, aes(x = NAME_2, y = elevation)) +
geom_boxplot() +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45))
#> Warning: Removed 51 rows containing non-finite outside the scale range
#> (`stat_boxplot()`).
# or just use terra::rasterize() + terra::boxplot():
y <- rasterize(v, r, "NAME_2")
boxplot(r,y)
Created on 2024-08-27 with reprex v2.1.1
Чрезвычайно хорошее предложение. Мне удалось воспроизвести это, используя файлы примеров terra. Однако я получил ошибку при использовании собственного файла формы и растра. > # извлечь, включить имена из векторных атрибутов > e <- extract(r,v) Ошибка в h(simpleError(msg, call)) : Fehler bei der Auswertung des Argumentes 'x' bei der Methodenauswahl für Funktion 'addAttrToGeom': coords должна быть матрицей из двух столбцов > extract(r, v) Ошибка в h(simpleError(msg, call)) : Fehler bei der Auswertung des Argumentes 'x' bei der Methodenauswahl für Funktion 'addAttrToGeom': координаты должны быть двух- матрица столбцов >
Я вижу, что мой шейп-файл отличается от образца Terra. поэтому e<-extract(r,v) или rasterize() не удалось > v Простая коллекция объектов с 12 объектами и 14 полями Тип геометрии: MULTIPOLYGON Размер: XY, XYZ Ограничивающая рамка: xmin: 2485410 ymin: 1075268 xmax: 2833842 ymax : 1295934 z_range: zmin: 0 zmax: 0 Прогнозируемый CRS: CH1903+ / LV95 # Тиббл: 12 × 15 RegionNumm RegionName Unterregio Unterreg_1 ObjNummer Version Shape_Leng Shape_Area DERegionNa FRRegionNa ITRegionNa DEBioBedeu <dbl> <chr> <dbl> <chr> хр> <дата> <двухместный> <дбл> –
Боюсь, вам придется предоставить полностью воспроизводимый пример вместе с наборами данных для менее общего ответа. Судя по вашему последнему комментарию, у вас также есть объект sf от sf::st_read(), а не SpatVector, который вы могли бы получить от terra::vect(). Под «terra или sf + stars» я имел в виду использование terra ИЛИ использование sf с starts. Но если речь идет о мультиполигональных фигурах, вы можете разделить их на многоугольники с помощью terra::disagg() .
Как мне правильно читать шейп-файл, а не read_sf? существует ли векторная функция? > v <- read_sf("......N2020_Revision_BiogeoRegion.shp") > str(v) sf [12 × 15] (S3: sf/tbl_df/tbl/data.frame) $ RegionNumm: num [1:12 ] 1 2 2 2 2 3 3 4 5 5 ... $ Unterregio: число [1:12] 11 21 22 23 24 31 32 41 51 52 .. $ Shape_Leng: число [1:12] 725117 334364 .. .. $ : num [1:23632, 1:2] 2641684 264 ..$ :Список 1 ..- attr(, "класс")= chr [1:3] "XY" "MULTIPOLYGON" "sfg" - attr( , "sf_column")= chr "geometry" - attr(*, "agr")= Коэффициент с 3 уровнями "константа", "агрегат",
как и в моем ответе, v <- vect("yourshapefilepath.shp")
Дальнейший вопрос. Что вы подразумеваете под идентификатором. e$NAME_2 <- v$NAME_2[e$ID] В моем шейп-файле у меня есть имена: имена: RegionNumm RegionName Unterregio Unterreg_1 ObjNummer Версия Shape_Leng Shape_Area DERegionNa FRRegionNa (и еще 4) type: <int> <chr> <int> < chr> <chr> <chr> <num> <num> <chr> <chr> Вы выбрали NAME_2, я выберу «Unterregio». В вашем файле я вижу ID_1, ID_2 как дополнительные имена. Я не уверен, какой идентификатор в вашем примере?
Столбец ID в результирующем кадре добавлен extract() для обозначения используемого многоугольника. Подробности можно проверить ?terra:::extract. В остальном, пожалуйста, отредактируйте свой вопрос так, чтобы проблемы, с которыми вы сталкиваетесь, можно было воспроизвести для других, начиная с загрузки библиотек и фактического минимального кода, который воспроизводит сбой на нашей стороне (мой ответ не загружается raster, но у вас есть RasterLayer). И, конечно же, нам также понадобятся наборы данных для запуска вашего кода. Подробнее о MRE — stackoverflow.com/collectives/r-language/articles/76995406/… .
Большое спасибо. Да, я полностью понимаю совет. Я все еще пытаюсь создать MRE на основе моего конкретного шейп-файла. Уже существует MRE на основе файлов примеров terra, и я не вижу никакой разницы по сравнению с моим шейп-файлом. До сих пор я думаю, что помимо MRE у меня есть координаты z, поэтому при загрузке я получил предупреждение. Вероятно, правильная загрузка шейп-файла устранит ошибку Extract(). # Загрузить файл shp > v <- vect("C:/Users/…._BiogeoRegion.shp") Предупреждающее сообщение: [vect] Координаты Z игнорируются > str(v) S4 class 'SpatVector' [пакет "terra"]
Решение прочитайте растровый файл с помощью terra::rast() и векторный шейп-файл с помощью vect(), затем вы можете следовать хорошему коду margusl
Чтобы указать, что ответ вам помог, вы можете пометить его как «принятый». Тенденция принимать/не принимать также видна через ваш профиль SO - stackoverflow.com/users/26852042/sibylle-st%c3%b6ckli - где соотношение вопросов/принятых вопросов может повлиять на мотивацию и выбор людей, которые могли бы потенциально ответить на ваши будущие вопросы. Можно отвечать (и принимать) на свои вопросы, но если контент в разделе ответов не совсем ответ, а скорее комментарий, он, вероятно, будет удален - stackoverflow.com/help/deleted-answers .
Проверьте свой
dataобъект. Не только имена/столбцы, но и реальный контент. Обратите внимание, что из включенного вами кода не совсем понятно, что именно было выполнено или вывод консоли был позже отредактирован вручную. Например, было ли фактическое выражениеextract(raster_file, shape_file)(без присвоения) илиextracted_values <- extract(raster_file, shape_file)