Criando arquivos WKT (Well Know Text)

Compreendendo arquivos WKT e suas aplicações

ogr / gdal / python / geoprocessing

Motivação

Tudo começou quando estava tentando abrir um arquivo .shp via GeoPandas (pacote Python) em um McBook (OS BigSur). Ao fazer a leitura do arquivo o seguinte erro era retornado.

...
IllegalArgumentException: geometries must not contain null elements
Shell is not a LinearRing

Nota: estava testando um McBook, uso sistema operacional Debian, e até o momento não tive este problema.

Após testar outros pacotes, transformar o arquivo para outros formatos, criar outros ambientes Python, entre outras tentativas, o erro persistia. Como não podia mais perder tempo, e o objetivo era inserir a bendita poligonal dentro do GeoDataframe, resolvi partir para o arquivo .wkt.

Introdução

Normalmente ao pensar em arquivos para objetos vetoriais (pontos, polígonos, etc), as primeiras opções que surgem na mente são: shp,GeoJSON, GeoPackage, etc. Porém, em alguns casos, usar o arquivo .wkt pode deixar tudo mais prático e rápido.

A extensão de arquivo .Wkt significa texto bem conhecido, e refere-se à um padrão de arquivo do OGC (Open Geospatial Consortium). Usualmente é utilizada para representar geometrias espaciais como texto em sistemas geoespaciais. De forma simplista, o arquivo .wkt contém uma representação de texto (string) para representar objetos geométricos vetoriais em um mapa. Informações extras estão disponíveis aqui.

Como já mencionado, com os arquivos wkt é possível representar objetos vetoriais, segue alguns exemplos:

Objeto WKT (string)
Pontos “POINT(3 2)”
Multi pontos “MULTIPOINT(3 3,5 8)”
Linha “LINESTRING(2.5 3.4,4.3 5.4)”
Polígono “POLYGON((0 0,1 4,3 4,3 2,0 0))”

Estudo de caso

Bibliotecas utilizados:

# ------------------------------------------------------------
# packages
import ogr
import geojson
import geopandas as gpd
import matplotlib.pyplot as plt

Para ler o arquivo .shp fiz uso da biblioteca ogr:

# ------------------------------------------------------------
# reading shape

shapePath = './field_011.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()
    
>>> ------------------------------------------------------- 
./field_011.shp has been opened

Uma vez que o arquivo foi aberto com sucesso, e o objeto layer recebeu a camada do arquivo field_011.shp, um novo objeto denominado feat recebeu a feição 0 do objeto layer. Por fim, este último foi “fechado” (recebeu None):

# get feature
feat = layer.GetFeature(0)
# close layer
layer = None

A próxima etapa foi obter as coordenadas do polígono, e criar o arquivo .wkt. O modo mais prático seria usar a biblioteca shapely, mas como a referida biblioteca também gerava o mesmo erro, optei por criar um objeto GeoJson com a poligonal em questão, e extrair as coordenadas para criar “manualmente” um .wkt.

  • Criando e lendo o .geoJson:
# create and read GeoJson
featGeoJson = feat.ExportToJson()
layerGeoJson = geojson.loads(featGeoJson)
>>> print(layerGeoJson)
{"geometry": {"coordinates": [[[-57.669078, -13.909359], [-57.670381, -13.912523], [-57.671333, -13.91486], [-57.672349, -13.917289], [-57.672533, -13.91768], [-57.672744, -13.917664], 
...
[-57.670604, -13.9091], [-57.669345, -13.909243], [-57.669109, -13.909274], [-57.669078, -13.909359]]], "type": "Polygon"}, "id": 0, "properties": {"Area_ha": 127.01, "Field": "011", "id": 1}, "type": "Feature"}
  • Extraindo as coordenadas do polígono e armazenando os valores em Tuple:
po = layerGeoJson["geometry"]["coordinates"]
coords = tuple([(x[0], x[1]) for x in po[0]])

>>> print(coords)
((-57.669078, -13.909359), (-57.670381, -13.912523), (-57.671333, -13.91486), (-57.672349, -13.917289), (-57.672533, -13.91768), (-57.672744, -13.917664), 
...
(-57.670604, -13.9091), (-57.669345, -13.909243), (-57.669109, -13.909274), (-57.669078, -13.909359))
  • Gerando o wkt:
# write a WKT from tuple
wktPol = "POLYGON (("
for idx in range(len(coords)):
    a, b = coords[idx][0], coords[idx][1]
    if idx == len(coords) - 1:
        wktPol = wktPol + str(str(a) + " " + str(b))
    else:
        wktPol = wktPol + str(str(a) + " " + str(b)) + ", "
wktPol = wktPol + "))"

>>> print(type(wktPol))
<class 'str'>
>>> print(wktPol)
POLYGON ((-57.669078 -13.909359, -57.670381 -13.912523, -57.671333 -13.91486, -57.672349 -13.917289, -57.672533 -13.91768, -57.672744 -13.917664,
...
-57.670604 -13.9091, -57.669345 -13.909243, -57.669109 -13.909274, -57.669078 -13.909359))

Para que o GeoPandas leia o .wkt, passamos o mesmo para um objeto GeoSeries e posteriormente cria-se o GeodataFrame:

gs = gpd.GeoSeries.from_wkt([wktPol])
gdf = gpd.GeoDataFrame(geometry=gs, crs="EPSG:4326")
>>> print(gdf)
                                            geometry
0  POLYGON ((-57.66908 -13.90936, -57.67038 -13.9...

Se necessário, basta criar outra coluna para identificar o polígono:

>>> gdf['id'] = ['field_11']
>>> print(gdf)
                                            geometry        id
0  POLYGON ((-57.66908 -13.90936, -57.67038 -13.9...  field_11

Como complemento, para gerar o .wkt via shapely, basta importar a referida biblioteca, criar um objeto com geometria polygon com as coordenadas (coords), e usar o método .to_wkt():

import shapely

shp = shapely.geometry.Polygon(coords)
wkt = shp.to_wkt()

 Share!

 
comments powered by Disqus