Otomatik NDVI analizi
https://geojson.io/#map=6/39.147/32.805 sitesinden incelemek istedi?imiz yeri seçiyoruz
bu indirdi?imiz dosyay? kullanarak hangi tarihler aras?nda veri almak istiyorsak o tarihler aras?nda filtreleme sayesinde veri topluyoruz
Stajyer Batuhan Ökmen bu verileri indiriyoruz ve ar?ivden ç?kar?yoruzbelirlenen nir ve red band?na bölüyoruz NDVI analizi yaparak tekrar bir resim olu?turuyoruz
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 4 11:02:01 2020
sentinentel uydu veri indirme
"""
from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
from datetime import date
import rasterio
import numpy as np
import pandas as pd
import os
#G?R??
api = SentinelAPI('KULLANICI ADI', '??FRE.', 'https://scihub.copernicus.eu/dhus')
#geojson dosyas? ile hangi alan? incelemek istiyorsak o alana gidiyoruz ve seçiyoruz kaydete basarak geojson olarak almak istedi?imizi söylüyoruz
footprint = geojson_to_wkt(read_geojson('map6.geojson'))
products = api.query(footprint,
date=('20191219', date(2019, 12, 29)),
platformname='Sentinel-2')
# pandas dataframe yap
products_df = api.to_dataframe(products)
# filtreleme
products_df_sorted = products_df.sort_values(['cloudcoverpercentage', 'ingestiondate'], ascending=[True, True])
products_df_sorted = products_df_sorted.head(5)
# ar?iv indirmesini burada yap?yoruz
api.download_all(products_df_sorted.index)
#index verisinden kullan?c? kodu çekme
veri_cekme=products_df_sorted.index
veri_cekme1=veri_cekme[0]
veri_cekme2=veri_cekme[1]
#title çekme i?lemi yapt?k
"""
Bu i?lem ar?ivden ç?karmak için gerekli ar?ivin ad? indirdi?imiz verinin title ad? oluyor
"""
arsiv_adi=api.get_product_odata(veri_cekme1)
arsiv_adi=arsiv_adi["title"]
arsiv_adi=str(arsiv_adi)
#ar?ivden ç?karmak için ar?iv kütüphanesini ekledik
from archive import Archive
a = Archive(arsiv_adi+'.zip')
a.extract()
img_data_klasor_ismi=os.listdir((arsiv_adi+".SAFE"+'/GRANULE'))
img_data_klasor_ismi=img_data_klasor_ismi[0]
img_data_klasor_ismi=str(img_data_klasor_ismi)
dosya_yer_=(arsiv_adi+".SAFE"+'/GRANULE/'+img_data_klasor_ismi+'/IMG_DATA')
resim_isim=os.listdir(dosya_yer_)
if resim_isim == "R10m" or "R20m" or "R60m":
dosya_yer_=(arsiv_adi+".SAFE"+'/GRANULE/'+img_data_klasor_ismi+'/IMG_DATA/R60m')
resim_isim=os.listdir(dosya_yer_)
resim_isim[2]
resim_isim[3]
jp2ler = [resim_isim[2],resim_isim[3]]
bands = []
#buras? bizim jp2 dosyalar?m?z? okuyacak
for jp2 in jp2ler:
with rasterio.open(dosya_yer_+"/"+jp2) as f:
bands.append(f.read(1))
#resimlerimizi ayr??t?rd?k özel bantlara
band_red=bands[0]
band_nir=bands[1]
else:
resim_isim[2]
resim_isim[3]
jp2ler = [resim_isim[2],resim_isim[3]]
bands = []
#buras? bizim jp2 dosyalar?m?z? okuyacak
for jp2 in jp2ler:
with rasterio.open(dosya_yer_+"/"+jp2) as f:
bands.append(f.read(1))
#resimlerimizi ayr??t?rd?k özel bantlara
band_red=bands[0]
band_nir=bands[1]
# Klasik NDVI denklemi ile hesaplama
np.seterr(divide='ignore', invalid='ignore')
# Calculate NDVI. This is the equation at the top of this guide expressed in code
ndvi = (band_nir.astype(float) - band_red.astype(float)) / (band_nir + band_red)
#su için yap?yorum bu analizi
##
###
ndvi=(band_red.astype(float) - band_nir.astype(float)) / (band_red + band_nir)
###
###
np.nanmin(ndvi), np.nanmax(ndvi)
# görüntümüze bakal?m renklerine ayr?lm?? bir görüntümüz var
# çizim yapaca??z bunun için gerekli kütüphaneler ekleniyor
import matplotlib.pyplot as plt
import matplotlib.colors as colors
# NDVI bilindi?i üzere 1 ve -1 aras?ndaki de?erlerde s?n?fland?r?l?r.
# Biz de bu de?erleri renklerle göstermek istiyoruz.
# Bunun için al?nan say?sal de?erleri farkl? renk spektrumlar?na atayarak elimizde NDVI için renklendirilmi? bir görüntümüz olacakt?r
#
# Bir orta nokta belirledik ve bu sola ve sa?a olacak ?ekilde renklendiriyoru renk spekturumunu da a?a??da payla?aca??m
class RenkNormalizasyonu(colors.Normalize):
def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
self.midpoint = midpoint
colors.Normalize.__init__(self, vmin, vmax, clip)
def __call__(self, value, clip=None):
x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
return np.ma.masked_array(np.interp(value, x, y), np.isnan(value))
min=np.nanmin(ndvi)
max=np.nanmax(ndvi)
mid=0.1
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
cmap = plt.cm.RdYlGn
cax = ax.imshow(ndvi, cmap=cmap, clim=(min, max), norm=RenkNormalizasyonu(midpoint=mid,vmin=min, vmax=max))
ax.axis('off')
ax.set_title('NDVI görüntüsü', fontsize=18, fontweight='bold')
cbar = fig.colorbar(cax, orientation='horizontal', shrink=0.65)
fig_kaydet="ndvi_sonuc/"+resim_isim[2]+".tif"
fig.savefig(fig_kaydet, dpi=200, bbox_inches='tight', pad_inches=0.7)
plt.show()
fig2 = plt.figure(figsize=(10,10))
ax = fig2.add_subplot(111)
plt.title("NDVI Histogram", fontsize=18, fontweight='bold')
plt.xlabel("NDVI values", fontsize=14)
plt.ylabel("# pixels", fontsize=14)
x = ndvi[~np.isnan(ndvi)]
numBins = 20
ax.hist(x,numBins,color='black',alpha=1)
histog_kaydet="histogram_sonuc/"+resim_isim[2]+".png"
fig2.savefig(histog_kaydet, dpi=200, bbox_inches='tight', pad_inches=0.7)
plt.show()
"""
2. part burada
burada görüntü i?leme ile ne durumda oldu?unu bulaca??z
"""
import cv2
#buras? 1 kere yap?lmas? gerekiyor
resim = cv2.imread(histog_kaydet,111)
cv2.imwrite("skysattif/donusumler//goruntu_isleme/resim_yenimis.jpg",resim)
resim_yeni="skysattif/donusumler//goruntu_isleme/resim_yenimis.jpg"
import cv2
import numpy as np
resim = cv2.imread(resim_yeni)
#bölge seçmeyi yap?yorum
#buras? 1.seviye kurak alan oluyor
kesilmis_karek1= resim[0:157,0:161]
# buras? 2. seviye kurak alan
kesilmis_karek2= resim[0:238,0:161]
# buras? 3. seviye kurak alan
kesilmis_karek3= resim[0:422,0:161]
#buras? 1.seviye az ye?il alan oluyor
kesilmis_kareaz1= resim[0:157,162:277]
# buras? 2. seviye az ye?il alan
kesilmis_kareaz2= resim[0:238,162:277]
# buras? 3. seviye az ye?il alan
kesilmis_kareaz3= resim[0:422,162:277]
#buras? 1.seviye çok ye?il alan oluyor
kesilmis_karey1= resim[0:157,278:462]
# buras? 2. seviye az ye?il alan
kesilmis_karey2= resim[0:238,278:462]
# buras? 3. seviye az ye?il alan
kesilmis_karey3= resim[0:422,278:462]
#listemelemek için
listem_ort=[]
kesilmis_karek1_ort=np.mean(kesilmis_karek1)
listem_ort.append(kesilmis_karek1_ort)
kesilmis_karek2_ort=np.mean(kesilmis_karek2)
listem_ort.append(kesilmis_karek2_ort)
kesilmis_karek3_ort=np.mean(kesilmis_karek3)
listem_ort.append(kesilmis_karek3_ort)
kesilmis_kareaz1_ort=np.mean(kesilmis_kareaz1)
listem_ort.append(kesilmis_kareaz1_ort)
kesilmis_kareaz2_ort=np.mean(kesilmis_kareaz2)
listem_ort.append(kesilmis_kareaz2_ort)
kesilmis_kareaz3_ort=np.mean(kesilmis_kareaz3)
listem_ort.append(kesilmis_kareaz3_ort)
kesilmis_karey1_ort=np.mean(kesilmis_karey1)
listem_ort.append(kesilmis_karey1_ort)
kesilmis_karey2_ort=np.mean(kesilmis_karey2)
listem_ort.append(kesilmis_karey2_ort)
kesilmis_karey3_ort=np.mean(kesilmis_karey3)
listem_ort.append(kesilmis_karey3_ort)
# yo?unluk hesapland?
listem_ort.sort()
en_kucuk=listem_ort[0]
if en_kucuk==kesilmis_karek1_ort:
toprak_turu="kurakl?k seviyesi 1"
elif en_kucuk==kesilmis_karek2_ort:
toprak_turu="kurakl?k seviyesi 2"
elif en_kucuk==kesilmis_karek3_ort:
toprak_turu="kurakl?k seviyesi 3"
elif en_kucuk==kesilmis_kareaz1_ort:
toprak_turu="Az ye?il 1"
elif en_kucuk==kesilmis_kareaz2_ort:
toprak_turu="Az ye?il 2"
elif en_kucuk==kesilmis_kareaz3_ort:
toprak_turu="Az ye?il 3"
elif en_kucuk==kesilmis_karey1_ort:
toprak_turu="Ye?il 1"
elif en_kucuk==kesilmis_karey2_ort:
toprak_turu="Ye?il 2"
elif en_kucuk==kesilmis_karey3_ort:
toprak_turu="Ye?il 3"
print("1 en yüksek de?eri ifade eder ye?il 1 en fazla ye?illik, kurakl?k 1 en kurak alan anlam?na gelir")
print("yo?unluk olarak: %s"%toprak_turu)
#ortalama hesaplan?yor
#burda en yak?n kom?uyu buluyoruz
from heapq import nsmallest
ortalama_yogunluk=nsmallest(1, listem_ort, key=lambda x: abs(x-(np.sum(listem_ort)/9)))
if ortalama_yogunluk==kesilmis_karek1_ort:
yogunluk_tipi="kurakl?k seviyesi 1"
elif ortalama_yogunluk==kesilmis_karek2_ort:
yogunluk_tipi="kurakl?k seviyesi 2"
elif ortalama_yogunluk==kesilmis_karek3_ort:
yogunluk_tipi="kurakl?k seviyesi 3"
elif ortalama_yogunluk==kesilmis_kareaz1_ort:
yogunluk_tipi="Az ye?il 1"
elif ortalama_yogunluk==kesilmis_kareaz2_ort:
yogunluk_tipi="Az ye?il 2"
elif ortalama_yogunluk==kesilmis_kareaz3_ort:
yogunluk_tipi="Az ye?il 3"
elif ortalama_yogunluk==kesilmis_karey1_ort:
yogunluk_tipi="Ye?il 1"
elif ortalama_yogunluk==kesilmis_karey2_ort:
yogunluk_tipi="Ye?il 2"
elif ortalama_yogunluk==kesilmis_karey3_ort:
yogunluk_tipi="Ye?il 3"
print("1 en yüksek de?eri ifade eder ye?il 1 en fazla ye?illik, kurakl?k 1 en kurak alan anlam?na gelir")
print("ortalama olarak: %s"%yogunluk_tipi)
print("yo?unluk olarak: %s"%toprak_turu)
"""
Buras? bölünen alanlar? görmek için var yok say?labilir
"""
"""
#gösterme kesilmi? kare kurak alan
cv2.imshow("kesilmis_karek1",kesilmis_karek1)
cv2.imshow("kesilmis_karek2",kesilmis_karek2)
cv2.imshow("kesilmis_karek3",kesilmis_karek3)
#gösterme kesilmi? az ye?il alan
cv2.imshow("kesilmis_kareaz1",kesilmis_kareaz1)
cv2.imshow("kesilmis_kareaz2",kesilmis_kareaz2)
cv2.imshow("kesilmis_kareaz3",kesilmis_kareaz3)
#gösterme kesilmi? çok ye?il alan
cv2.imshow("kesilmis_karey1",kesilmis_karey1)
cv2.imshow("kesilmis_karey2",kesilmis_karey2)
cv2.imshow("kesilmis_karey3",kesilmis_karey3)
cv2.waitKey(0)
cv2. destroyAllWindows()
type(resim)
"""
# pasta grafi?ine dönü?türüyoruz
import matplotlib.pyplot as plt
%matplotlib inline
labels="kurakl?k seviyesi 1","kurakl?k seviyesi 2","kurakl?k seviyesi 3","Az ye?il 1","Az ye?il 2","Az ye?il 3","Ye?il 1","Ye?il 2","Ye?il 3"
sizes=[kesilmis_karek1_ort,kesilmis_karek2_ort,kesilmis_karek3_ort,kesilmis_kareaz1_ort,kesilmis_kareaz2_ort,kesilmis_kareaz3_ort,kesilmis_karey1_ort,kesilmis_karey2_ort,kesilmis_karey3_ort]
explode = (0,0,0,0.2,0.2,0.2,0,0.1,0)
figl, axl = plt.subplots()
axl.pie(sizes,explode=explode,labels=labels,autopct="%1.1f%%", shadow=True, startangle=90)
axl.axis("equal")
figl.savefig("pasta_grafik.png")
plt.show()
