# Segundo o método do valor informativo, produz um cartograma que representa a suscepetibilidade relativa à ocorrência de movimentos de vertente, em cada unidade territorial (célula) - GDAL
import os, shutil, numpy, math
from osgeo import gdal, ogr, osr
def split_gdal(entrada, w, classe, sre):
    # Abrir input layer
    inDriver = ogr.GetDriverByName("ESRI SHAPEFILE")
    inDataSource = inDriver.Open(entrada, 0)
    inLayer = inDataSource.GetLayer()
    inLayer.SetAttributeFilter("value = " + str(classe))
    # Criar saida
    outDriver = ogr.GetDriverByName("ESRI SHAPEFILE")
    outDataSource = outDriver.CreateDataSource(w + "\\parte_" + str(classe) + ".shp")
    outLayer_name = "parte_" + str(classe)
    outLayer = outDataSource.CreateLayer(outLayer_name, geom_type=ogr.wkbPolygon)
    # Criar campos
    inLayerDefn = inLayer.GetLayerDefn()
    for i in range(0, inLayerDefn.GetFieldCount()):
        fieldDefn = inLayerDefn.GetFieldDefn(i)
        fieldName = fieldDefn.GetName()
        outLayer.CreateField(fieldDefn)
    outLayerDefn = outLayer.GetLayerDefn()
    for inFeature in inLayer:
        outFeature = ogr.Feature(outLayerDefn)
        for i in range(0, outLayerDefn.GetFieldCount()):
            fieldDefn = outLayerDefn.GetFieldDefn(i)
            fieldName = fieldDefn.GetName()
            outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
        geom = inFeature.GetGeometryRef()
        outFeature.SetGeometry(geom.Clone())
        outLayer.CreateFeature(outFeature)
    inDataSource.Destroy()
    outDataSource.Destroy()
    cria_prj(sre, w + "\\parte_" + str(classe) + ".prj")
def intersection(point, polygon, out):
    drv = ogr.GetDriverByName("ESRI SHAPEFILE")
    point = drv.Open(point, 0)
    point_lyr = point.GetLayer()
    polygon = drv.Open(polygon, 0)
    polygon_lyr = polygon.GetLayer()
    union1 = ogr.Geometry(3)
    for feature in polygon_lyr:
        geom = feature.GetGeometryRef()
        union1 = union1.Union(geom)
    union2 = ogr.Geometry(point_lyr.GetGeomType())
    for feat in point_lyr:
        geom = feat.GetGeometryRef()
        union2 = union2.Union(geom)
    intersect = union1.Intersection(union2)
    dstshp = drv.CreateDataSource(out)
    dstlayer = dstshp.CreateLayer(out, geom_type=ogr.wkbPolygon)
    featureDefn = dstlayer.GetLayerDefn()
    i_poligonos = {}
    id_feature = 0
    for polygon in intersect:
        figura = str(polygon)
        l_polygon = str.split(figura[10:-2], ",")
        print l_polygon
        i_poligonos.update({id_feature:l_polygon})
        id_feature +=1
    for chave in i_poligonos.keys():
        feature = ogr.Feature(featureDefn)
        ring = ogr.Geometry(ogr.wkbLinearRing)
        for coord in i_poligonos[chave]:
            coord = str.split(coord, " ")
            try:
                x = float(coord[0])
                y = float(coord[1])
            except:
                x = float(coord[1][1:])
                y = float(coord[2])
            ring.AddPoint(x, y)
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        feature.SetGeometry(poly)
        dstlayer.CreateFeature(feature)
def pre_valor_informativo(dic, fenomeno, factor, workspace, resolucao, sre):
    c = 0
    d = {}
    for chave in dic.keys():
        intersection(fenomeno, dic[chave], workspace + "\\mov_" + str(factor) + "_" + str(c) + ".shp")
        cria_prj(sre, workspace + "\\mov_" + str(factor) + "_" + str(c) + ".prj")
        shp2rst(workspace + "\\mov_" + str(factor) + "_" + str(c) + ".shp", resolucao, -9999, workspace + "\\r_mov_" + str(factor) + "_" + str(c) + ".tif")
        try:
            count_pixeis_fenomeno = conta_ocorrencias(workspace + "\\r_mov_" + str(factor) + "_" + str(c) + ".tif")
        except:
            count_pixeis_fenomeno = 0
        d.update({int(chave):count_pixeis_fenomeno})
        c += 1
    return d
def reclass_vector(raster, dicionario, output, resolucao, sre):
    img = gdal.Open(raster)
    v = numpy.array(img.ReadAsArray())
    new = []
    lst_classes = []
    for chave in dicionario.keys():
        lst_classes.append(chave)
    for i in v:
        linha = []
        for e in i:
            for c in lst_classes:
                if int(e) == int(c):
                    cat = c
                else:
                    continue
            if e < 0 or e > 20:
                continue
            else:
                linha.append(dicionario[cat][1])
        new.append(linha)
    array2raster(raster, output, int(resolucao), int(resolucao), numpy.array(new[::-1]), sre)
def main_vi_GDAL(fenomeno, wfactores, output, sresp):
    for chave in d_factores.keys():
        # Converter o raster para vectorial
        rst2Polygon(d_factores[chave][0], w + "\\v_factor_" + str(chave) + ".shp")
        cria_prj(prj, w + "\\v_factor_" + str(chave) + ".prj")
        # Partir o factor em partes
        os.mkdir("C:\\s_factor" + str(chave))
        id_cls = id_classes(d_factores[chave][1])
        for classe in id_cls:
            split_gdal(w + "\\v_factor_" + str(chave) + ".shp", "C:\\s_factor" + str(chave), classe, prj)
        # Associar o id da classe ao ficheiro respectivo
        lst_files = lista_shapefile("C:\\s_factor" + str(chave))
        dic_cls = ass_cls_split(id_cls, lst_files)
        # Calcular o valor informativo de cada classe
        pre_vi = pre_valor_informativo(dic_cls, fenomeno, chave, w, cellsize, prj)
        _vi_ = {}
        px_classe = nr_pixel_classe(d_factores[chave][1])
        for ch in px_classe.keys():
            id_classe = int(ch)
            count_classe = px_classe[ch]
            count_ocorrencia = pre_vi[id_classe]
            val_info = (float(count_ocorrencia) / count_classe) / (float(area_ocorrencias) / area_total)
            if val_info == 0:
                lnvi = -2 * 1000000
            else:
                lnvi = ((math.log10(val_info)) * 1000000)
            _vi_.update({id_classe:[val_info, lnvi]})
        # Reclassificar o raster original, atribuindo-lhe os valores do valor informativo
        print _vi_
        reclass_vector(d_factores[chave][0], _vi_, w + "\\recls_factor_" + str(chave) + ".tif", cellsize, sresp)
        lst_reclss.append(w + "\\recls_factor_" + str(chave) + ".tif")
    # somar todos os raster
    vectores = []
    for i in lst_reclss:
        vector = create_array(i)
        vectores.append(vector)
    final = []
    for i in range(len(vectores[0])):
        linha = []
        for e in range(len(vectores[0][i])):
            l_vi = []
            for u in range(len(vectores)):
                vi = vectores[u][i][e]
                l_vi.append(vi)
            _vi_ = 0
            for vi in l_vi:
                _vi_ += vi
            linha.append(_vi_)
        final.append(linha)