R ile Topografik Veri İndirme 2 - DEM Verisi İndirme ve Görselleştirme

R
gis
r-spatial
base R
R ile Topografik Topografik Veri İndirme 2
Author

Mehmet Göktuğ Öztürk

Published

March 10, 2025

Merhabalar, blogun yeni 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.

Geçtiğimiz hafta, üzerinden uzunca bir zaman geçtikten sonra, blogun üçüncü yazısını yayımlamıştım. geodata paketini kullanarak SRTM verisi indirmiş ve görselleştirmiştim. SRTM, orijinal olarak 30 metre çözünürlüklü bir veri ancak geodata paketi bu veriyi, maksimum 30 arc second (yaklaşık 1 km) çözünürlüklü olarak sağlıyor. Bu çözünürlük, çoğu görselleştirme ve ekolojik modelleme çalışması için oldukça yeterli, özellikle de Türkiye gibi büyük coğrafyalarda çalıştığımız zaman. Ancak çalışmamızın konusu, amacı ya da çalışma alanının boyutuna göre, topografyayı daha iyi temsil edecek, daha yüksek çözünürlüklü verilere ihtiyaç duyabiliriz. Geçtiğimiz yazının devamı niteliğindeki bu yazıda, R kullanarak yüksek çözünürlüklü SRTM verilerine nasıl erişebileceğimizi anlatacağım. Yüksek çözünürlü veri indireceğimiz için Türkiye ölçeğinde veri indirmek ve o veriyi işlemek daha fazla zaman alacaktır. Bu sebeple bu yazı İzmir ölçeğinde olacak.

Yazının Akışı

  1. Paketlerin yüklenmesi
  2. DEM verisinin indirilmesi
  3. Verinin görselleştirilmesi

1. Paketlerin yüklenmesi

Bu yazı için 4 paket kullanacağız. Mekânsal vektör veriler için sf, raster veriler için terra, mülki idare sınırları için yine rgeoboundaries ve interaktif görselleştirme için de mapview paketini kullanacağız. Bu yazıdaki yeni paketimiz elevatr. DEM verilerine erişmek için oldukça pratik bir paket.

library(sf)
library(terra)
library(elevatr)
library(rgeoboundaries)
library(mapview)

2. DEM verisinin indirilmesi

elevatr paketini kullanarak birden fazla kaynaktan (SRTM, ALOS 3D, Amazon) veri indirmek mümkün. Yukarıda da bahsettiğim gibi bu yazıda SRTM verisi indireceğiz. Bunun için ilk olarak Türkiye mülki idare sınırı verisini indirelim.

tr_il <- gb_adm1("TUR")
tr_il
Simple feature collection with 81 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 25.66545 ymin: 35.80768 xmax: 44.81766 ymax: 42.1048
Geodetic CRS:  WGS 84
First 10 features:
   shapeGroup shapeType      shapeName shapeISO                 shapeID
1         TUR      ADM1          Adana    TR-01 25984515B95172477822815
2         TUR      ADM1       Adıyaman    TR-02 25984515B90872828599679
3         TUR      ADM1 Afyonkarahisar    TR-03 25984515B26209284550223
4         TUR      ADM1           Ağrı    TR-04 25984515B39003465173278
5         TUR      ADM1         Amasya    TR-05 25984515B32583172380009
6         TUR      ADM1        Antalya    TR-07 25984515B97476213604692
7         TUR      ADM1         Artvin    TR-08  25984515B1602304411378
8         TUR      ADM1          Aydın    TR-09 25984515B63739470532168
9         TUR      ADM1      Balıkesir    TR-10 25984515B55384476443375
10        TUR      ADM1         Ankara    TR-06 25984515B47806653651907
   shapeCanonical                       geometry
1        province MULTIPOLYGON (((35.38791 36...
2        province MULTIPOLYGON (((37.861 37.4...
3        province MULTIPOLYGON (((30.48061 38...
4        province MULTIPOLYGON (((43.77542 39...
5        province MULTIPOLYGON (((36.3878 40....
6        province MULTIPOLYGON (((30.42781 36...
7        province MULTIPOLYGON (((41.87016 40...
8        province MULTIPOLYGON (((27.33197 37...
9        province MULTIPOLYGON (((26.7224 39....
10       province MULTIPOLYGON (((31.96673 38...
plot(st_geometry(tr_il))

Şekil 1: Türkiye mülki idare sınırları.

Şimdi de İzmir sınırlarını çekelim.

izm <- tr_il[tr_il$shapeName == "İzmir", ]
izm
Simple feature collection with 1 feature and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 26.19417 ymin: 37.81523 xmax: 28.49308 ymax: 39.38549
Geodetic CRS:  WGS 84
   shapeGroup shapeType shapeName shapeISO                 shapeID
41        TUR      ADM1     İzmir    TR-35 25984515B59477790330217
   shapeCanonical                       geometry
41       province MULTIPOLYGON (((28.48507 38...
plot(st_geometry(izm))

Şekil 2: İzmir mülki idare sınırları.

Bir de interaktif olarak kontrol edelim.

mapview(izm)

Şekil 3: İzmir mülki idare sınırları.

İzmir sınırları hazır olduğuna göre DEM verisini indirebiliriz.

elevatr paketi, “aws” hariç diğer kaynaklardan veri indirmek için OpenTopography API’sini kullanıyor. Bunun için siteye üye olmalı ve API Key oluşturmalısınız. Bu, işlem ücretsiz ve oldukça basit. Sonrasında set_opentopo_key() fonksiyonunu kullanarak key’inizi girmelisiniz. Ardından R oturumunu yeniden başlattığınızda verilere erişebilirsiniz.

İzmir sınırlarını kullanarak DEM verisini indirmek için get_elev_raster() fonksiyonunu kullanacağız. İzmir DEM verisini indirmek için locations = izm, yüksek çözünürlüklü SRTM verisi için src = "gl1" ve veriyi İzmir sınırlarına göre kırpmak için de clip = "locations" argümanını kullanacağız.

dem <- get_elev_raster(izm, src = "gl1", clip = "locations")
dem
class      : RasterLayer 
dimensions : 5653, 8276, 46784228  (nrow, ncol, ncell)
resolution : 0.0002777778, 0.0002777778  (x, y)
extent     : 26.19431, 28.49319, 37.81514, 39.38542  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : file90761e26a847 
values     : -25, 2148  (min, max)
plot(dem)

Şekil 4: İzmir topografya haritası, raster paketinin varsayılan plot() fonksiyonu.

Yüksek çözünürlüklü DEM verisini indirdik. Gördüğünüz gibi veri RasterLayer sınıfında. Bunun sebebi elevatr paketinin hâlen raster paketini kullanması. Veriyi terra’nın SpatRaster sınıfına çeviriyorum. Ardından da bir önceki yazıdaki gibi reclassify edip, 0’dan küçük olan tüm değerleri 0’a eşitleyeceğim. Bunu terra paketi içindeki classify() fonksiyonu ile yapabiliriz.

dem <- rast(dem)
dem <- classify(
  dem, 
  rcl = matrix(c(-Inf, 0, 0), ncol = 3, byrow = TRUE)
)
dem
class       : SpatRaster 
dimensions  : 5653, 8276, 1  (nrow, ncol, nlyr)
resolution  : 0.0002777778, 0.0002777778  (x, y)
extent      : 26.19431, 28.49319, 37.81514, 39.38542  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
source(s)   : memory
name        : file90761e26a847 
min value   :                0 
max value   :             2148 
plot(dem)

Şekil 5: İzmir topografya haritası, terra paketinin varsayılan plot() fonksiyonu.

Şimdi de görselleştirmeye geçebiliriz.

3. Verinin görselleştirilmesi

Veriyi görselleştirmek için yine base R kullanacağım. Bu sefer daha iyi bir renk paleti seçiyor ve kuzey oku ile ölçek ekliyorum.

breaks <- seq(0, 2150, 10)                                                    # kirilimlari onceden belirliyoruz

par(family = "Montserrat", cex = 1.4)                                         # font ailesi ve boyutu
plot(
  dem,                                                                        # raster veri 
  breaks = breaks,                                                            # kirilimlar
  type = "continuous",                                                        # kirilimlari manuel belirledigimiz icin lejanti manuel olarak cont. yapmamiz gerekiyor
  col = tidyterra::hypso.colors(                                              # tidyterra paketiyle palet seciyorum  
    length(breaks) - 1,                                                       # renk sayisinin kırilim boyutundan bir az olmasi icin
    palette = "wiki-schwarzwald-cont"
  ), 
  las = 1,                                                                    # y eksenindeki yaziların yatay olmasi icin
  maxcell = prod(dim(dem)),                                                   # maksimum cozunurlukte cizmek icin
  plg = list(at = seq(0, 2100, 250), tic = "out"),                            # lejant ayarlari
  mar = c(0.5, 1.5, 0, 4)                                                     # marjin 
)
plot(st_geometry(izm), add = TRUE)                                            # izmir sinirlari
north(c(26.35, 39.2), type = 1, label = expression(bold("KUZEY")), lwd = 1.2) # kuzey oku
sbar(10, "bottomleft", below = "km", adj = c(0.5, -1), lonlat = TRUE)         # olcek
mtext("Veri: SRTM", side = 1, adj = 0.99, line = 1.7)                         # kaynak eklemek icin

Şekil 6: İzmir topografya haritası, terra ve base plot ile topografya haritası.

Geçenki yazıdaki histogramı şimdi de İzmir için çizelim, rakımın İzmir’deki dağılımına bir bakalım.

breaks <- seq(0, 2150, 50)

par(
  mar = c(5, 5.2, 1, 1),
  mgp = c(3.6, 0.7, 0),                                                     # etiket uzakligini ayarliyoruz
  las = 1, 
  family = "Montserrat", 
  bty = "l",                                                                # box'ın tipini L seklinde yapmak icin
  cex.lab = 1.2,
  xpd = NA                                                                  # etiketi grafik sinirlari disina tasimaya izin verir
)
hist(dem, xlab = "Rakım [m]", ylab = "Sıklık", main = "", breaks = breaks)  # histogram
axis(side = 1, at = breaks, tck = -0.030 / 4, labels = FALSE)               # x eksenindeki kucuk tikler
box()                                                                       # eksenlere kutu cizmek icin

Şekil 7: İzmir’de rakımın dağılımı.

Grafik, sergilediği sağa çarpık dağılımla İzmir’in bir kıyı kenti olduğunu oldukça net bir şekilde gösteriyor.

Bir sonraki yazıda görüşürüz.

Bilimle ve huzurla kalınız!