library(raster)
library(terra)
library(sf)

################ Формируем папку с необходимыми данными и присваиваем короткие имена########
Data = "/home/taras/GRIDS"
setwd(Data)
list.files(pattern="\\.TIF$")

B1 = 0.0000275*(raster ("LT05_L2SP_180028_19940304_20200913_02_T2_SR_B1.TIF"))- 0.2
B2 = 0.0000275*(raster ("LT05_L2SP_180028_19940304_20200913_02_T2_SR_B2.TIF"))- 0.2
B3 = 0.0000275*(raster ("LT05_L2SP_180028_19940304_20200913_02_T2_SR_B3.TIF"))- 0.2
B4 = 0.0000275*(raster ("LT05_L2SP_180028_19940304_20200913_02_T2_SR_B4.TIF"))- 0.2
B5 = 0.0000275*(raster ("LT05_L2SP_180028_19940304_20200913_02_T2_SR_B5.TIF"))- 0.2
TEMP_C = (0.00341802 * (raster ("LT05_L2SP_180028_19940304_20200913_02_T2_ST_B6.TIF"))+ 149.0)-273.15
B7 = 0.0000275*(raster ("LT05_L2SP_180028_19940304_20200913_02_T2_SR_B7.TIF"))- 0.2
CLOUD = raster ("CLOUD_04_94.TIF")
#Проверяем систему координат снимков
B1 
################Читаем вектор################

dnestr_area <- st_read("/home/taras/GRIDS/Dnestr_area/Dnestr_area.shp")

#plot(dnestr_area, main = "Dnestr area", axes = TRUE, border = "blue")

################ВЕГЕТАТИВНЫЕ ИНДЕКСЫ########################
                    #NDVI
#NDVI = (NIR − Red)/(NIR + Red)

NDVI_C = (B4-B3)/(B4+B3)
NDVI = NDVI_C*CLOUD
#plot(NDVI)
#plot(dnestr_area, col="transparent", add=T)
NDVI_crop = crop(NDVI, dnestr_area)
mask_NDVI_crop = mask(NDVI_crop,  dnestr_area)

writeRaster(mask_NDVI_crop, "/home/taras/GRIDS/NDVI_03_94.tiff")

rm("NDVI_C", "NDVI", "NDVI_crop", "mask_NDVI_crop")

NDVI_max <- raster("~/GRIDS/NDVI_03_94.tiff")
#N_min <- NDVI_max@data@min
N_max <- NDVI_max@data@max

                    #VCI
#VCI=(NDVI - NDVI_min) / (NDVI_max - NDVI_min)
VCI=(NDVI_max - 0.2) / (N_max - 0.2)
writeRaster(VCI, "/home/taras/GRIDS/VCI_03_94.tiff")


                    #TEMPERATURA  Landsat 5 and Landsat 7 Level-2
#TEMP_C=(0.00341802 * B6 + 149.0)-273.15# Считаю сразу
TEMP = TEMP_C*CLOUD
TEMP_crop = crop(TEMP, dnestr_area)
mask_TEMP_crop = mask(TEMP_crop, dnestr_area)
writeRaster(mask_TEMP_crop, "/home/taras/GRIDS/TEMP_03_94.tiff")

rm("TEMP_C", "TEMP", "TEMP_crop", "mask_TEMP_crop")

TEMP_max <- raster("~/GRIDS/TEMP_03_94.tiff")
T_min <- TEMP_max@data@min
T_max <- TEMP_max@data@max

                    #TCI
#TCI = (BTmax - BT) / (BTmax – BTmin) 
TCI = (T_max - TEMP_max) / (T_max - T_min)
writeRaster(TCI, "/home/taras/GRIDS/TCI_03_94.tiff")

                    #VHI
#VHI = (VCI+TCI)/2
VHI = (VCI+TCI)/2
writeRaster(VHI, "/home/taras/GRIDS/VHI_03_94.tiff")

rm("TEMP_max", "TCI", "VHI", "VCI")

                    #Green Chlorophyll Index (CI green)
#GCI = (NIR / Green) - 1

GCI_C = (B4 / B2) - 1
GCI = GCI_C*CLOUD
GCI_crop = crop(GCI, dnestr_area)
mask_GCI_crop = mask(GCI_crop,  dnestr_area)
writeRaster(mask_GCI_crop, "/home/taras/GRIDS/GCI_03_94.tiff")

rm("GCI_C", "GCI", "GCI_crop", "mask_GCI_crop")

                    #SLAVI
#SLAVI = NIR/(RED+SWIR)

SLAVI_C = B4/(B3+B7)
SLAVI = SLAVI_C*CLOUD
SLAVI_crop = crop(SLAVI, dnestr_area)
mask_SLAVI_crop = mask(SLAVI_crop,  dnestr_area)
writeRaster(mask_SLAVI_crop, "/home/taras/GRIDS/SLAVI_03_94.tiff")

rm("SLAVI_C", "SLAVI", "SLAVI_crop", "mask_SLAVI_crop")

                    #SIPI
#SIPI = (NIR - BLUE) / (NIR - RED)

SIPI_C = (B4 - B1) / (B4 - B3)
SIPI = SIPI_C*CLOUD
SIPI_crop = crop(SIPI, dnestr_area)
mask_SIPI_crop = mask(SIPI_crop,  dnestr_area)
writeRaster(mask_SIPI_crop, "/home/taras/GRIDS/SIPI_03_94.tiff")

rm("SIPI_C", "SIPI", "SIPI_crop", "mask_SIPI_crop")

################ВОДНЫЕ ИНДЕКСЫ########################
               #MSI
#MSI = SWIRI/NIR

MSI_C = B5/B4
MSI = MSI_C*CLOUD
MSI_crop = crop(MSI, dnestr_area)
mask_MSI_crop = mask(MSI_crop,  dnestr_area)
writeRaster(mask_MSI_crop, "/home/taras/GRIDS/MSI_03_94.tiff")

rm("MSI_C", "MSI", "MSI_crop", "mask_MSI_crop")


               #NDMI
#NDMI = (NIR - SWIR1) / (NIR + SWIR1)

NDMI_C = (B4 - B5) / (B4 + B5)
NDMI = NDMI_C*CLOUD
NDMI_crop = crop(NDMI, dnestr_area)
mask_NDMI_crop = mask(NDMI_crop,  dnestr_area)
writeRaster(mask_NDMI_crop, "/home/taras/GRIDS/NDMI_03_94.tiff")

rm("NDMI_C", "NDMI", "NDMI_crop", "mask_NDMI_crop", "CLOUD")
rm("B1", "B2", "B3", "B4", "B5", "B7")
gc()