R • Ekoloji • Veri Bilimi
10 Ara 2023
Merhabalar, blogun ikinci yazısına hoş geldiniz. Takıldığınız ve anlamadığınız yerler olursa lütfen yorum yapmaya çekinmeyiniz. Ayrıca katkılarınızı ve eleştirilerinizi de bekliyorum. Keyifli okumalar.
Blogun ilk yazısını yazarken, ikincisi de birkaç haftaya çıkar, nasıl olsa kodlar hazır diye düşünüyordum ancak tezdir, iştir derken uzadı. Yoğunluk yetmiyormuş gibi bilgisayarın format ihtiyacı ve NTFS dosya sistemine sahip harici diskin linux ile uyumsuzluğu da eklenince yazının taslağı ve kodlarını kaybettim. Sonuç olarak -birbirinden güzel bahanelerimi özetleyecek olursam- gecikmiş, kaybedilmiş ve sıfırdan yazılmış bir yazıyla karşınızdayım.
Bu blog yazısı, bir öncekinin devamı niteliğinde sayılabilir. R ile eBird verilerini işlemeye devam edeceğiz. İlk yazıda birinci soru üzerine odaklanmıştık, bunda ise ikinci soruya odaklanacağız:
Temel R ve GIS bilgisine sahip olmanız yazıdaki birçok yeri anlamanız için yeterli olacaktır. Ancak bu alanlarda yeniyseniz de endişelenmeyin. Anlamadığınız her şeyi yorum olarak sorabilirsiniz. Blogun da ana amaçlarından birisi bu zaten. Ayrıca kuşlara ilginiz olmasa bile, bu yazının bazı temel mekânsal analizleri öğrenmek için faydalı olacağını düşünüyorum.
| Paket | Açıklama |
|---|---|
| tidyverse | Çoğunlukla veri manipülasyonu ve görselleştirme üzerine paketler içeren bir paket koleksiyonu |
| sf | Simple Features: Mekânsal vektör verileri işlemek için |
| terra | Raster verileri işlemek için |
| tidyterra | terra paketi için tidyverse metotları sağlayan bir paket |
| mapview | İnteraktif haritalar için |
| rgeoboundaries | Mülki idare sınırlarını indirmek için |
Bu paketler kurulu değilse aşağıdaki kod bloğu ile kurabilirsiniz. Bu kod bloğu, paketi R’a yüklemeye çalışacak, eğer yükleyemezse kuracaktır. Paketler kuruluysa, library() fonksiyonu ile yükleyebiliriz.
# eger paketler yuklu degilse onlari kur
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("sf")) install.packages("sf")
if (!require("terra")) install.packages("terra")
if (!require("tidyterra")) install.packages("tidyterra")
if (!require("mapview")) install.packages("mapview")
if (!require("rgeoboundaries")) remotes::install_github("wmgeolab/rgeoboundaries")
# paketleri R'a yukle
library(tidyverse) # bircok veri isini kolaylastirmak için
library(sf) # r'da mekansal vektor verileri islemek icin
library(terra) # r'da mekansal raster verileri islemek icin
library(tidyterra) # terra ile tidyverse metotlari kullanmak icin
library(mapview) # interaktif haritalar cizmek icin
library(rgeoboundaries) # tr ilce sinirlarina erismek icinİkinci soruya cevap verebilmek için kuş, ilçe sınırları ve dijital yükseklik modeline (DEM) ihtiyacımız var. İlk ilçeden başlayalım. Kuş ve DEM verisini işlemek için ilçeye ihtiyaç duyacağız ne de olsa.
İlk yazıda olduğu gibi bu sefer de geoBoundaries verisini kullanacağız. Şimdi ilçe verilerini indirip, içinden Kızılcahamam’ı çıkartalım.
tr_ilce <- gb_adm2(country = "Turkey", type = "SSCGS") # type = "SSCGS" argumaniyla basitlestirilmiş versiyonunu indiriyoruz
khamam <- tr_ilce[tr_ilce$shapeName == "Kızılcahamam", ]
khamamSimple feature collection with 1 feature and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 32.27539 ymin: 40.21093 xmax: 32.91882 ymax: 40.75389
Geodetic CRS: WGS 84
shapeGroup shapeType shapeName shapeISO shapeID
388 TUR ADM2 Kızılcahamam <NA> 54988432B2925959455693
shapeCanonical geometry
388 Districts MULTIPOLYGON (((32.51627 40...Kızılcahamam’ı doğru bir şekilde aldık gibi görünüyor. Verimizi {mapview} paketini kullanarak görselleştirelim, bir kontrol edelim. Bu paket, interaktif haritalar yapmak için oldukça işlevsel.
mapview(khamam)Gördüğünüz gibi her şey yolunda. :)
Kara akbaba - Aegypius monachus, Avrasya’da geniş dağılıma sahip bir yırtıcı türü. 3 metreye varan kanat açıklığıyla en büyük kuşlardan birisi ve azalan bir popülasyona sahip. Türkiye, türün en önemli yaşam alanlarından birisi. Tür, özellikle Kuzey Anadolu’da, orman - bozkır geçişindeki yaşlı ormanları tercih ediyor. Bu alanlardan birisi de bu yazıda incelediğimiz yer olan Kızılcahamam. Kızılcahamam’ın yaşlı karaçam - Pinus nigra ormanları, türün üremesi için oldukça elverişli.
Türle alakalı kısa bir bilgi verdikten sonra verilere tekrardan dönebiliriz. Kuş verisini nasıl elde ettiğimize ilk yazıda değinmiştik. Tekrardan uzun uzun anlatmaya gerek yok.
Öncelikle kuş verimizi R’a yükleyelim:
ebird <- read_delim("./ebird/ebd_TR_relApr-2023.txt")print(ebird)# A tibble: 2,403,720 × 50
`GLOBAL UNIQUE IDENTIFIER` `LAST EDITED DATE` `TAXONOMIC ORDER` CATEGORY
<chr> <dttm> <dbl> <chr>
1 URN:CornellLabOfOrnithology:E… 2021-04-15 12:59:44 10043 species
2 URN:CornellLabOfOrnithology:E… 2021-04-15 12:58:56 5625 species
3 URN:CornellLabOfOrnithology:E… 2021-04-15 13:04:32 1361 species
4 URN:CornellLabOfOrnithology:E… 2021-04-14 00:03:13 22215 species
5 URN:CornellLabOfOrnithology:E… 2021-04-05 11:01:32 5787 species
6 URN:CornellLabOfOrnithology:E… 2021-04-21 02:26:13 29239 species
7 URN:CornellLabOfOrnithology:E… 2021-04-14 00:03:35 22156 species
8 URN:CornellLabOfOrnithology:E… 2021-04-05 11:01:32 291 species
9 URN:CornellLabOfOrnithology:E… 2018-09-20 02:46:53 5945 species
10 URN:CornellLabOfOrnithology:E… 2018-09-20 02:46:53 5945 species
# ℹ 2,403,710 more rows
# ℹ 46 more variables: `TAXON CONCEPT ID` <chr>, `COMMON NAME` <chr>,
# `SCIENTIFIC NAME` <chr>, `SUBSPECIES COMMON NAME` <chr>,
# `SUBSPECIES SCIENTIFIC NAME` <chr>, `EXOTIC CODE` <chr>,
# `OBSERVATION COUNT` <chr>, `BREEDING CODE` <chr>,
# `BREEDING CATEGORY` <chr>, `BEHAVIOR CODE` <chr>, `AGE/SEX` <chr>,
# COUNTRY <chr>, `COUNTRY CODE` <chr>, STATE <chr>, `STATE CODE` <chr>, …Gördüğünüz gibi 2,403,720 gözlem (satır) ve 50 değişkene (sütun) sahip bir veri tablomuz var. Bu verilerin çoğuna ihtiyaç duymuyoruz. Dolayısıyla işimize yaramayanlardan kurtulmamız işimizi çokça kolaylaştıracaktır. Sadece işimize yarayan sütunları seçip, kara akbaba türüne ait gözlemleri içerecek şekilde filtreleyelim.
karaakbaba <- ebird |>
select(6, 7, 11, 12, 29, 30) |> # burada indeks kullanarak sectik, sutun isimleriyle de secebiliriz
filter(`SCIENTIFIC NAME` == "Aegypius monachus")
karaakbaba# A tibble: 1,910 × 6
`COMMON NAME` `SCIENTIFIC NAME` `OBSERVATION COUNT` `BREEDING CODE` LATITUDE
<chr> <chr> <chr> <chr> <dbl>
1 Cinereous Vul… Aegypius monachus 1 <NA> 40.5
2 Cinereous Vul… Aegypius monachus 5 <NA> 40.5
3 Cinereous Vul… Aegypius monachus 1 <NA> 40.0
4 Cinereous Vul… Aegypius monachus 3 <NA> 40.5
5 Cinereous Vul… Aegypius monachus 1 <NA> 40.5
6 Cinereous Vul… Aegypius monachus 1 <NA> 40.5
7 Cinereous Vul… Aegypius monachus 2 <NA> 40.5
8 Cinereous Vul… Aegypius monachus 2 <NA> 40.5
9 Cinereous Vul… Aegypius monachus 1 <NA> 40.5
10 Cinereous Vul… Aegypius monachus 2 <NA> 40.5
# ℹ 1,900 more rows
# ℹ 1 more variable: LONGITUDE <dbl>Şimdi değişken sayımız 6’ya, gözlem sayımız ise 1910’a düştü. Veri, çok daha yönetilebilir bir hâlde. Bu etkiyi, değişken isimlerini değiştirerek daha da artırabiliriz. Daha kısa, tamamı küçük harflerden oluşan, latinize ve boşluklar yerine alt çizgi kullandığımız isimler daha kullanışlı olacaktır. Herhangi bir hatayla karşılaşma riskimiz düşecektir. Bunun için ilk olarak karaakbaba data frame’inin değişken (sütun) isimlerine bakalım. names() ya da colnames() fonksiyonunu kullanabiliriz.
names(karaakbaba)[1] "COMMON NAME" "SCIENTIFIC NAME" "OBSERVATION COUNT"
[4] "BREEDING CODE" "LATITUDE" "LONGITUDE" Şimdi de yeni isimlerden oluşan bir karakter vektörü oluşturalım ve bunu karaakbaba data frame’inin değişken isimlerine atayalım.
names(karaakbaba) <- c(
"eng_name", "sci_name", "obs_count", "breeding_code", "y", "x"
)
names(karaakbaba)[1] "eng_name" "sci_name" "obs_count" "breeding_code"
[5] "y" "x" Gördüğünüz gibi değişken isimlerini istediğimiz gibi değiştirdik.
Artık ilk yüklediğimiz eBird verisini R’dan silebiliriz. Bu veri, 2 milyondan fazla gözlemden oluştuğu için RAM’in şişmesine ve R’ın çökmesine yol açabilir. Ardından da çöp toplayıcıyı (garbage collector) çalıştırmak faydalı olacaktır. R, bu işi genelde otomatik yapıyor ancak büyük veriler silindiği zaman gc()’nin çalıştırılması öneriliyor.
rm(ebird)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 3964934 211.8 9938292 530.8 7650458 408.6
Vcells 8520332 65.1 125684295 958.9 145568976 1110.7Şimdi de mekânsallaştırmaya geçebiliriz. Verimiz şu anda bir data frame. Dolayısıyla onunla mekânsal işler yapmak için verimizi mekânsallaştırmamız gerekiyor. Bunun için {sf} paketini kullanacağız. Bu paket, mekânsal vektör veriler işlemek için geliştirmiş, OGC (Open Geospatial Consortium) standardı bir paket. Ben de işlerimin büyük bir kısmında onu tercih ediyorum.
Koordinatların olduğu sütunları ve koordinat sistemini tanımlayarak verimizi sf objesine dönüştürelim.
karaakbaba_sf <- st_as_sf(
karaakbaba, coords = c("x", "y"), crs = "EPSG:4326"
)
karaakbaba_sfSimple feature collection with 1910 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 26.50652 ymin: 36.02635 xmax: 44.5648 ymax: 42.0362
Geodetic CRS: WGS 84
# A tibble: 1,910 × 5
eng_name sci_name obs_count breeding_code geometry
* <chr> <chr> <chr> <chr> <POINT [°]>
1 Cinereous Vulture Aegypius… 1 <NA> (32.62995 40.46758)
2 Cinereous Vulture Aegypius… 5 <NA> (32.62995 40.46758)
3 Cinereous Vulture Aegypius… 1 <NA> (34.61423 40.01914)
4 Cinereous Vulture Aegypius… 3 <NA> (32.62995 40.46758)
5 Cinereous Vulture Aegypius… 1 <NA> (32.6622 40.48013)
6 Cinereous Vulture Aegypius… 1 <NA> (32.62403 40.4557)
7 Cinereous Vulture Aegypius… 2 <NA> (32.62403 40.4557)
8 Cinereous Vulture Aegypius… 2 <NA> (32.62403 40.4557)
9 Cinereous Vulture Aegypius… 1 <NA> (32.62403 40.4557)
10 Cinereous Vulture Aegypius… 2 <NA> (32.62403 40.4557)
# ℹ 1,900 more rowsKara akbaba verisi artık mekânsal bir formatta. Geometri tipi, koordinat sistemi gibi öznitelikler barındırıyor.
Veriyi daha iyi tanımak için ilk görselleştirmeleri yapalım. Acaba kara akbabalar Türkiye’de nasıl dağılıyor? Bunun için tekrardan mapview() fonksiyonunu kullanıyoruz.
mapview(karaakbaba_sf)Bu işlem de tamam. Gördüğünüz gibi ülkemizin kuzeyinde yoğunlaşan bir popülasyon mevcut.
Sıra geldi kuş verimizin Kızılcahamam’a göre alt kümesini almaya.
karaakbaba_subset <- karaakbaba_sf[khamam, ]
karaakbaba_subsetSimple feature collection with 139 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 32.29039 ymin: 40.24768 xmax: 32.8661 ymax: 40.64056
Geodetic CRS: WGS 84
# A tibble: 139 × 5
eng_name sci_name obs_count breeding_code geometry
<chr> <chr> <chr> <chr> <POINT [°]>
1 Cinereous Vulture Aegypius… 1 <NA> (32.62995 40.46758)
2 Cinereous Vulture Aegypius… 5 <NA> (32.62995 40.46758)
3 Cinereous Vulture Aegypius… 3 <NA> (32.62995 40.46758)
4 Cinereous Vulture Aegypius… 1 <NA> (32.6622 40.48013)
5 Cinereous Vulture Aegypius… 1 <NA> (32.62403 40.4557)
6 Cinereous Vulture Aegypius… 2 <NA> (32.62403 40.4557)
7 Cinereous Vulture Aegypius… 2 <NA> (32.62403 40.4557)
8 Cinereous Vulture Aegypius… 1 <NA> (32.62403 40.4557)
9 Cinereous Vulture Aegypius… 2 <NA> (32.62403 40.4557)
10 Cinereous Vulture Aegypius… 2 <NA> (32.62403 40.4557)
# ℹ 129 more rowsKızılcahamam’da, eBird veritabanında kayıtlı 139 tane kara akbaba kaydı varmış. Bunlara bir göz atıp DEM verisiyle ilgilenmeye başlayalım. Kara akbaba dağılımıyla alakalı asıl haritayı en son yapacağız.
plot(
st_geometry(khamam),
main = "Kızılcahamam'da Kara Akbaba Dağılımı",
lwd = 2,
border = "grey30",
axes = TRUE,
reset = FALSE
)
plot(
st_geometry(karaakbaba_subset),
cex = 1,
pch = 3,
lwd = 1.5,
add = TRUE
)
Hemen hızlıca mapview ile de bakalım.
mapview(khamam, lwd = 3, alpha.regions = 0) + mapview(karaakbaba_subset)Kara akbabanın Kızılcahamam’daki dağılımını ortaya çıkarttık. Şimdi de bazı çevresel değişkenlerle ilişkisi nasılmış, -kabaca- ona bakalım. Ben, topografik değişkenlerle olan ilişkisini merak ettim ve bu sebeple DEM verisine bakacağız.
DEM (Digital Elevation Model), kısaca, yeryüzünün yükseklik bilgisini sayısal olarak temsil eden modele verilen isim. Her bir grid için, sahip olduğu yükseklik verisini barındırır. Topografik ve hidrolojik birçok analizde sıklıkla kullanılıyor. Rakımın biyolojik sistemlerdeki önemi nedeniyle, ekoloji ve biyocoğrafyada da yaygın bir kullanımı var.
Burada COP-DEM (Copernicus DEM) kullandım. Veriye şu bağlantıdan ulaşabilirsiniz. Belki daha sonra bu verinin nereden, nasıl indirilebileceğiyle alakalı bir başka yazı yazarım.
DEM verisi, raster formatında bir veri. Raster veriyi, bir tür resim ya da bir tür matris olarak düşünebilirsiniz. Bu veri, hücrelerden - gridlerden oluşur ve her bir grid için bir değer içerir. Genellikle, yükseklik, eğim, sıcaklık, yağış gibi sürekli veriler bu formattadır.
Biz bu veri için {terra} paketini kullanacağız. Oldukça hızlı ve işlevsel bir paket.
# dem verisini yukleyelim
dem <- rast("./dem_4326.tif")
demclass : SpatRaster
size : 7556, 22983, 1 (nrow, ncol, nlyr)
resolution : 0.0008333333, 0.0008333333 (x, y)
extent : 25.66542, 44.81792, 35.80792, 42.10458 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : dem_4326.tif
name : Copernicus_DSM_COG_30_N35_00_E025_00_DEM
min value : -86.01817
max value : 5107.41357 Veriyi yazdırdığımızda, çözünürlüğünden, dört köşesinin koordinatına ve koordinat sistemine kadar birçok detayı görebiliyoruz.
plot(dem)
Kolayca haritasını da çizebiliyoruz. Gördüğünüz gibi Türkiye topografik haritası.
Şimdi verimizi Kızılcahamam ilçe sınırlarına göre maskeleyelim. khamam, sf formatında olduğu için vect() ile {terra}’nın kendi vektör veri tipi olan SpatVector’e dönüştürüyoruz.
dem_subset <- crop(dem, vect(khamam), mask = TRUE)
dem_subsetclass : SpatRaster
size : 651, 772, 1 (nrow, ncol, nlyr)
resolution : 0.0008333333, 0.0008333333 (x, y)
extent : 32.27542, 32.91875, 40.21125, 40.75375 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s) : memory
varname : dem_4326
name : Copernicus_DSM_COG_30_N35_00_E025_00_DEM
min value : 749.5305
max value : 2069.5696 plot(dem_subset)
terrain() ile eğim ve bakıyı hesaplayalım.
slp <- terrain(dem_subset, "slope")
asp <- terrain(dem_subset, "aspect")par(mfrow = c(1, 2))
plot(slp, main = "Kızılcahamam Eğim Haritası")
plot(asp, main = "Kızılcahamam Bakı Haritası")
Şimdi de her bir kara akbaba gözlem noktasının sahip olduğu çevresel değişkenleri çıkartalım ve ardından karaakbaba verisiyle birleştirelim.
karaakbaba_subset$ID <- 1:nrow(karaakbaba_subset)
ext_dem <- extract(dem, vect(karaakbaba_subset))
ext_slp <- extract(slp, vect(karaakbaba_subset))
ext_asp <- extract(asp, vect(karaakbaba_subset))
ext_dem <- left_join(karaakbaba_subset, ext_dem, by = "ID")
ext_slp <- left_join(ext_dem, ext_slp, by = "ID")
ext_asp <- left_join(ext_slp, ext_asp, by = "ID")
karaakbaba_ext <- ext_asp
names(karaakbaba_ext)[7] <- "altitude"
glimpse(karaakbaba_ext)Rows: 139
Columns: 9
$ eng_name <chr> "Cinereous Vulture", "Cinereous Vulture", "Cinereous Vul…
$ sci_name <chr> "Aegypius monachus", "Aegypius monachus", "Aegypius mona…
$ obs_count <chr> "1", "5", "3", "1", "1", "2", "2", "1", "2", "2", "2", "…
$ breeding_code <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ geometry <POINT [°]> POINT (32.62995 40.46758), POINT (32.62995 40.4675…
$ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
$ altitude <dbl> 1058.8070, 1058.8070, 1058.8070, 1045.1802, 1136.2458, 1…
$ slope <dbl> 11.188035, 11.188035, 11.188035, 22.408733, 11.157285, 1…
$ aspect <dbl> 309.599090, 309.599090, 309.599090, 242.915925, 7.856463…Gördüğünüz gibi kara akbaba verimizin içinde artık her bir nokta için rakım, eğim ve bakı bilgisi mevcut. Bu aşamadan sonra modellerle çalışabiliriz ancak bu, bu yazının konusu değil. Ayrıca elimizdeki veri bunun için ne kadar yeterli, ayrı bir tartışma konusu. Şimdi çıkarttığımız topografik verilere bir göz atıp asıl haritamızı yapmaya geçebiliriz.
karaakbaba_ext |>
st_drop_geometry() |>
pivot_longer(cols = c("altitude", "slope", "aspect")) |>
ggplot(aes(x = value, fill = name, colour = "name")) +
geom_histogram(alpha = .6) +
scale_fill_viridis_d() +
scale_color_viridis_d() +
theme_minimal() +
theme(
legend.position="none",
) +
facet_wrap(~name, scales = "free_x")
Uzun bir yazı oldu yine ancak sonuna geldik. Şimdi Kızılcahamam’daki kara akbaba gözlemlerinin dağılımını, DEM ile birlikte görselleştirmeye çalışalım.
ggplot() + # grafigi baslatiyor
geom_spatraster( # dem verisini ekliyoruz
data = dem_subset,
na.rm = TRUE
) +
geom_sf( # kizilcahamam ilce sinirlarini ekliyoruz
data = khamam,
fill = NA
) +
geom_sf( # kara akbaba gozlem verisini ekliyoruz
data = karaakbaba_subset,
aes()
) +
scale_fill_wiki_c() + # dem verisi icin renk paleti ekliyoruz
labs( # etiketleri yaziyoruz
title = "Kara Akbaba Dağılım Haritası",
subtitle = "Kızılcahamam'daki kara akbaba gözlemlerinin dağılımı",
fill = "Yükseklik (m)",
caption = "Veri: eBird | COP-DEM"
) +
theme_minimal() + # tema seciyoruz
theme( # temanin ozelliklerini belirliyoruz
plot.background = element_rect("white", colour = "white"),
text = element_text(family = "ubuntu mono")
)
Güzel oldu sanki. :)
Bilimle ve huzurla kalınız.