Módulo 01

Objetivo

Ler, editar e gravar arquivo vetorial (shape).

Material e métodos

Linux Debian Buster
GDAL 2.4.0, released 2018/12/14

Na manipulação e processamento de aquivo vetorial, utilizaremos um arquivo shape com os municípios do estado de MT (Figura 1).

Figura 1: Arquivo shape do estado de MT com Tabela de Atributos

Leitura e gravação de arquivo vetorial (shape)

Importação dos pacotes

Inicialmente importar os pacotes

from osgeo import ogr, osr, gdal

Leitura do arquivo

Para abrir o arquivo vetorial, define-se o caminho e o parâmetro 0 ou 1:

# ------------------------------------------------------------
# reading shape
# ------------------------------------------------------------
shapePath = './shape/mt_municipios_2020.shp'
ds = ogr.Open(shapePath, 1)  # 0 means read-only. 1 means writeable

# Check to see if shapefile was found.
if ds is None:
    print("------------------------------------------------------- \n")
    print('Could not open %s' % (shapePath))
    print('Check shapefile')
else:
    print("------------------------------------------------------- \n")
    print('%s' ' has been opened' % (shapePath))
    layer = ds.GetLayer()

Nesse bloco foi especificado o caminho que o arquivo encontra-se, e o método Open da classe OGR é executado com o argumento 1 para possibilitar a edição do arquivo. Na sequência, criou-se uma verificação da leitura do arquivo. Se o arquivo foi aberto com sucesso, o objeto layer recebe a classe Layer via ds.GetLayer(). Verifique aqui todos os métodos disponíveis na classe Layer aqui

Para conhecer o número de feições do arquivo e o CRS do aquirvo utiliza-se:

# Get the number of features in the layer
numFeatures = layer.GetFeatureCount()
print(numFeatures)
>>> print(numFeatures)
141

# get CRS
>>> print(layer.GetSpatialRef())
GEOGCS["SIRGAS 2000",
    DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6674"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4674"]]

No exemplo utilizado, há 141 feições, que no caso representam os municípios do estado. Para obter os campos (atributos ou nome das colunas) contidos na Tabela de Atributos há diversas opções, vamos listar duas:

# opt 01
# get name of atributes in shape file
>>> ldefn = layer.GetLayerDefn()
>>> ldefn.GetFieldCount()  # number of columns
5

lsAtribNames = []
for idx in range(ldefn.GetFieldCount()):
    lsAtribNames.append(ldefn.GetFieldDefn(idx).name)

>>> print(lsAtribNames)
['CD_MUN', 'NM_MUN', 'SIGLA_UF', 'AREA_KM2', 'AREA_ha']


# opt 02
lsAtribNames = [field.name for field in layer.schema]
>>> print(lsAtribNames)
['CD_MUN', 'NM_MUN', 'SIGLA_UF', 'AREA_KM2', 'AREA_ha']

Como resultado verifica-se que o arquivo tem 5 campos: “CD_MUN”, “NM_MUN”, “SIGLA_UF”, “AREA_KM2”, “AREA_ha”. Agora. caso deseje-se obter os valores dos campos de uma dada feição:

# ------------------------------------------------------------
# opt 01
# get feature
for feat in layer.GetFeature(0):
    print(feat)

>>> 5100102
Acorizal
MT
850.763
85076.296

for feat in layer.GetFeature(1):
    print(feat)
>>> 5100201
Água Boa
MT
7549.233
754928.252

# ------------------------------------------------------------
opt 02
>>> feat = layer.GetFeature(0)
>>> feat.GetField('NM_MUN')
'Acorizal'
>>> feat.GetField('SIGLA_UF')
'MT'
>>> feat.GetField('AREA_KM2')
850.763
>>> feat.GetField('AREA_ha')
85076.296

Note que ao mudar o número entre parênteses, obtém-se valores de outra feição. Agora vamos editar o arquivo shape, como exercício vamos criar um campo (atributo) para o cálculo da área em m².

Edição do arquivo

A lógica consiste em criar o novo campo desejado e posteriormente preencher com os valores desejados. No caso é á área em m².

# new field
newField = ogr.FieldDefn("Area_m2", ogr.OFTReal)
newField.SetWidth(16)
newField.SetPrecision(1)
layer.CreateField(newField)

# calculate areaM2
for idx in range(layer.GetFeatureCount()):
    feat = layer.GetFeature(idx)
    areaM2 = feat.GetField('AREA_ha') * 10000
    feat.SetField("Area_m2", areaM2)
    layer.SetFeature(feat)

# close to write file
layer = feat = None

A variável newField é um objeto da classe osgeo.ogr.FieldDefn e recebe o “comprimento” do número e a precisão das casas decimais. Na sequência, no objeto layer é criado o campo Area_m2. Note que é muito similar a lógica da calculadora de campo do QGIS. Posteriormente o referido campo é preenchido com os valores da área em m². Por fim, é necessário fechar as variáveis para o arquivo ser gravado em disco.

Se verificarmos a Tabela de atributo do arquivo no QGIS, nota-se a coluna com a área em m² (Figura 2). Não esqueça de importar o arquivo shape novamente.

The power of GDAL and OGR

Leitura, edição e gravação de arquivo vetorial (shape)