Python ile Sayısal Yükseklik Modeli Analizleri

Bu atölyede SRTM raster verileri işlenerek bazı temel analizler gerçekleştirilecektir.

Kullanılan kütüphaneler: rasterio, geopandas, rasterstats

• rasterio: mekansal raster veri işlemleri için
• geopandas: mekansal veri operasyonları için
• rasterstats: vektör geometriler ile raster verilerden istatistik çıkarmak için
• contextily: internetten altık haritaları çekmek için

1. Uygulama: SRTM Sayısal Yükseklik Verisi Okuma ve Bölgesel İstatistik

• https://search.earthdata.nasa.gov/search/granules?p=C1546314043-LPDAAC_ECS adresinden NASADEM Merged DEM Global 1 arc second V001 veri seti araması ile sağ taraftan indireceğimiz bölgenin alanını kare içine alarak arama yapıyoruz ve ilgili bölgenin .zip uzantılı .hgt yükseklik veri setini indiriyoruz.

• İndirdiğimiz bölgeye ait ilçe sınırlarını vektör veri formatında https://gadm.org/ sitesinden indiriyoruz. Level-2 olan veri ilçeleri içermektedir. GeoJSON formatında Level-2 veri bu çalışmada yeterli olacaktır.

• Level-2 Türkiye GeoJSON Verisi: https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_TUR_2.json.zip

• İndirdiğimiz ilçe verilerini okuyup görselleştirelim.

# read turkey level-2 geojson
ilceler = gpd.read_file('gadm41_TUR_2.json')
ilceler.plot()

• İstanbul, Bursa ve Kocaeli illerini seçelim

secileniller = ilceler[(ilceler['NAME_1'] == 'Istanbul') | (ilceler['NAME_1'] == 'Bursa') | (ilceler['NAME_1'] == 'Kocaeli')]

• Bursa'nın, Karamürsel, Orhangazi ve İznik ilçelerini seçelim.

bursailceler = ilceler[(ilceler['NAME_2'] == 'Karamürsel') | (ilceler['NAME_2'] == 'Orhangazi') | (ilceler['NAME_2'] == 'İznik')]

• Seçilen ilçeleri çizdirelim. Bu ilçelere göre sayısal yükseklik modelimizi keseceğiz.

bursailceler.plot()

• Şimdi sayısal yükseklik modeli verimizi öncelikle rasterio ile okuyup, rioshow ile gösterip, contextily'dan alacağımız arkaplan harita ile sunalım. Üzerine'de seçtiğimiz illerin katmanı nı ekleyelim.

r = rasterio.open("n40e029.hgt")

f, ax = plt.subplots(1, figsize=(9, 9))
ax = rioshow(r, alpha=0.5, zorder=2, ax=ax)
cx.add_basemap(ax, crs=r.crs) 
secileniller.plot(facecolor="none", edgecolor="#F93822", zorder=1, ax=ax)

• Seçtiğimiz Bursa ilçeleri için görüntüyü kesmeden önce bir kere daha görselleştirelim.

r = rasterio.open("n40e029.hgt")

f, ax = plt.subplots(1, figsize=(9, 9))
ax = rioshow(r, alpha=0.5, zorder=2, ax=ax)
cx.add_basemap(ax, crs=r.crs) 
bursailceler.plot(facecolor="none", edgecolor="#F93822", zorder=1, ax=ax)

• Şimdi bölgesel istatistik hesabı için rasterstats kütüphanesinden zonal_stats fonksiyonu ile ilçelerimiz için min, max, median gibi değerleri hesaplayalım.

from rasterstats import zonal_stats

elevations2 = zonal_stats(bursailceler, 'n40e029.hgt',  stats="count min mean max median")
elevations2 = pandas.DataFrame(elevations2)

• Son aşamada sayısal yükseklik verisini bu ilçe sınırlarına göre keserek bir bir .tif uzantılı dosya üreteceğiz.

• Bunun için öncelikle ilçeler verimizi geojson olarak kaydedelim.

bursailceler.to_file("bursa.geojson", driver="GeoJSON")

• rio mask komutu ile görüntümüzü bu vektör dosya sınırına göre keserek yeni bir dosya yazdıralım.

! rio mask n40e029.hgt \
           bursa.tif \
           --crop \
           --geojson-mask\
           bursa.geojson

• Sonrasında görselleştirelim.

f, ax = plt.subplots(1, figsize=(9, 9))
r = rasterio.open("bursa.tif")
rioshow(r, zorder=1, alpha=0.75, ax=ax)
cx.add_basemap(ax, alpha=0.5, crs=r.crs)