#-----Inicio do Script (rodar funcoes antes)-----#
#-----Funcoes se encontram no script 'funcoes_NDVI'-----#
rodar_bibliotecas(TRUE)
raiz = 'C:\\Users\\Yury\\Desktop\\InCeres\\Projeto_Fertilidade\\Sprint\\02\\'

#-----Ler arquivos de entrada-----#
shapefile_contorno = readOGR(paste0(raiz, 'fields_6_areas_validas.shp'))
#imagem_satelite = stack(paste0(raiz, 'imagem.tif'))
imagem_satelite_04 = stack(paste0(raiz, 'IMG_04.tif'))
imagem_satelite_05 = stack(paste0(raiz, 'IMG_05.tif'))

#-----Transformar projecao e recortar .TIF-----#
shapefile_contorno = transformar_em_coordenadas_utm(shapefile_contorno)
datum_correto = proj4string(shapefile_contorno)

numero_talhoes = 1:(length(shapefile_contorno[,1]))
#i=3
for (i in 1:length(numero_talhoes)){
  
  #-----Definir a separacao dos poligonos para clusterizar um de cada vez-----#
  shapefile_contorno_do_talhao_original = shapefile_contorno[i,1]
  shapefile_contorno_do_talhao = shapefile_contorno[i,1]
  shapefile_contorno_do_talhao = buffer(shapefile_contorno_do_talhao, width = 60)
  shapefile_contorno_para_recorte = reprojetar_para_recortar_utm(shapefile_contorno_do_talhao,imagem_satelite_04)
  imagem_recortada_04 = crop(imagem_satelite_04, shapefile_contorno_para_recorte)
  imagem_recortada_05 = crop(imagem_satelite_05, shapefile_contorno_para_recorte)
  reprojetado_04 = projectRaster(imagem_recortada_04, crs = datum_correto)
  reprojetado_05 = projectRaster(imagem_recortada_05, crs = datum_correto)
  
  
  #-----calcular NDVI-----#
  a = reprojetado_05+reprojetado_04
  b = reprojetado_05-reprojetado_04
  
  raster_ndvi = b/a
  raster_produtividade = 0.1665*exp(8.9245*(raster_ndvi))
  raster_extracao = raster_produtividade*0.63
  
  q2 = list(raster_ndvi, raster_produtividade, raster_extracao)
  tipo = c('NDVI', 'Produtividade', 'Extracao_K')
  #j=1
  for(j in 1:length(q2)){
    
    raster_ndvi_recortado = raster::mask(q2[[j]], shapefile_contorno_do_talhao)
    shapefile_pontos_ndvi = rasterToPoints(raster_ndvi_recortado, spatial = TRUE)
    shapefile_pontos_ndvi = intersect(shapefile_pontos_ndvi, shapefile_contorno_do_talhao_original)
    dataframe_ndvi = rasterToPoints(raster_ndvi_recortado, spatial = FALSE)
    write.table(dataframe_ndvi,(paste0(raiz, i, '_', tipo[j], '_dataframe_ndvi.txt')), sep = ';', row.names = F)
    
    #plot(shapefile_pontos_ndvi)
    
    #-----Carregando os dados para interpolacao-----#
    dados_ndvi = data.frame(dataframe_ndvi)
    
    attach(dados_ndvi)
    #names(dados_ndvi)
    nome_analise = 'Krigagem_NDVI'
    
    #-----Carrega dados no geoR (pacote de geoestatistica)-----#
    geo_dados_ndvi = read.geodata(paste0(raiz, i, '_', 'dataframe_ndvi.txt'), sep = ';', header = TRUE, coords.col=1:2, data.col=3)
    
    #-----Estatistica analitica-----#
    #estb = basicStats(layer, ci=0.95)
    #tks = ks.test(layer, 'pnorm', mean=mean(layer), sd=sd(layer))
    #tks
    #shapiro.test(layer)
    #jarqueberaTest(layer)
    #boxplot(layer)
    
    #-----Geoestatistica - escolhendo um melhor ajuste de modelo-----#
    #-----Distancia Max e Min-----#
    long_max = max(x)
    lat_max = max(y)
    mdta = sqrt(mx^2+my^2)
    mdta
    mdta/3
    
    #-----Gera os diferentes semivariogramas-----#
    #-----maxima distancia e numero de logs-----#
    mdt = 130
    nlg = 6 
    
    semi_variograma_teorico = variog(geo_dados_ndvi, max.dist = mdt, uvec = nlg)
    #semi_variograma_teorico1 = variog(geo_dados_ndvi, max.dist = mdt, uvec = nlg, direction = pi/8)
    #semi_variograma_teorico2 = variog(geo_dados_ndvi, max.dist = mdt, uvec = nlg, direction = pi/4)
    #semi_variograma_teorico3 = variog(geo_dados_ndvi, max.dist = mdt, uvec = nlg, direction = pi/2)
    
    h = semi_variograma_teorico$u
    v = semi_variograma_teorico$v
    npar = semi_variograma_teorico$n
    ltmx = (max(h)+0.4*max(h))
    ltmy = (max(v)+0.4*max(v))
    ic = (1.96*(sqrt(2*v)/sqrt(npar)))
    is = abs(v+ic)
    ii = abs(v-ic)
    
    #-----Ajuste dos semivariogramas teoricos-----#
    esferico = variofit(semi_variograma_teorico, cov.model='sph', max.dist=max(h), messages=F)
    #exponencial = variofit(semi_variograma_teorico, cov.model='exponential', max.dist=max(h), messages=F)
    #gaussiano = variofit(semi_variograma_teorico, cov.model='gaussian', max.dist=max(h),messages=F)
    
    #-----Separando os parametros dos modelos-----#
    #-----Esferico-----#
    c0.esf = esferico$nugget
    c.esf = esferico$cov.pars[1]
    a.esf = esferico$cov.pars[2]
    #-----Exponencial-----#
    #c0.exp = exponencial$nugget
    #c.exp = exponencial$cov.pars[1]
    #a.exp = exponencial$cov.pars[2]
    #-----Gaussiano-----#
    #c0.gau = gaussiano$nugget
    #c.gau = gaussiano$cov.pars[1]
    #a.gau = gaussiano$cov.pars[2]
    
    #-----Gera o grid para receber os valores interpolados-----#
    tamanho_do_pixel = 20 # Tamanho do pixel para interpolação
    range_x = as.integer(range(x))
    range_y = as.integer(range(y))
    grid.map = expand.grid(x=seq(from=range_x[1], to=range_x[2], by = (tamanho_do_pixel)),
                           y=seq(from=range_y[1], to=range_y[2], by = (tamanho_do_pixel)))
    
    #-----Krigagem-----#
    semi_variograma_ajustado = esferico
    krg_01 = krige.conv(geo_dados_ndvi, locations = grid.map, krige = krige.control(obj.model = semi_variograma_ajustado))
    length(grid.map[,1])
    fim = data.frame(grid.map$x, grid.map$y, krg_01$predict)
    coordenadas = data.frame(grid.map$x, grid.map$y)
    coordenadas = coordinates(coordenadas)
    dados = data.frame(fim$krg_01.predict)
    #write.table(fim, paste0(raiz, i, '_', 'teste_grid.txt'), sep = ',', row.names = F)
    shp_final = SpatialPointsDataFrame(coords = coordenadas, data = dados)
    shp_final = intersect(shp_final, shapefile_contorno_do_talhao_original)
    
    salvar = "C:\\Users\\Yury\\Desktop\\InCeres\\Projeto_Fertilidade\\Sprint\\02"
    writeOGR(shp_final,
             dsn = salvar,
             layer = paste0('talhao_',i, '_', tipo[j], '_NDVI_do_CACETE'),
             driver = 'ESRI Shapefile',
             overwrite_layer = TRUE)
    
    detach(dados_ndvi)
  }
}