Otimizando o Python: Cython

Otimizar um programa é fazer com que este programa ou algumas de suas funções rode mais rápido.Uma regra básica é que quanto menor o nível da linguagem de programação (quanto mais perto da linguagem de máquina) maior a possibilidade de otimização.

python

A linguagem Python é famosa por ser uma linguagem de alto nível e de alta produtividade por facilitar a vida do programador. A linguagem C é famosa por ser uma linguagem de baixo nível extremamente eficiente em termos computacionais.

Uma forma comum de otimizar o código Python fazer com que ele execute módulos escritos em C, a biblioteca padrão do Python vem com uma porção de módulos que funcionam assim, por exemplo o itertools. Algumas bibliotecas externas também usam este recurso para maior eficiência, o Numpy é um ótimo exemplo, ele é extremamente eficiente.

Para o usuário é possível escrever suas próprias extensões em C que podem ser chamadas diretamente pelo código Python, contudo isto requer o conhecimento avançado de C.

A biblioteca Cython elimina esta dificuldade, com ela é possivel escrever extensões em Python que são convertidas para uma forma mais eficiente em C, são compiladas e podem ser chamadas como módulos com o comando import.

Abaixo segue um exemplo onde a concatenação de duas strings é repetida 70 milhões de vezes, esta operação é realizada primeiro em Python e depois em um  módulo convertido para C, sobre como usar o Cython, consultar a documentação.

import cProfile
import pyximport
pyximport.install()
import Cythonteste

def juntar_texto(texto1, texto2):
    for i in range(70000000):
        texto_final = texto1+texto2

print "Normal: ----------------------------------------------------"
cProfile.run("juntar_texto('Feliz ', '2013')")
print "Cython: ----------------------------------------------------"
cProfile.run("Cythonteste.juntar_texto('Feliz ', '2013')")

O módulo Cythonteste consiste de um arquivo com exatamente a mesma função juntar_texto, ele é convertido e compilado automaticamente pelo módulo pyximport.
Para medir o tempo de execução é utilizado o cProfile. Veja o resultado, são gastos 8.51 segundos no Python contra 2.35 segundo no código compilado para rodar a mesma função.

Normal: ----------------------------------------------------
4 function calls in 8.509 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
 1 0.000 0.000 8.509 8.509 <string>:1(<module>)
 1 6.053 6.053 8.509 8.509 teste_comparativo.py:6(juntar_texto)
 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
 1 2.456 2.456 2.456 2.456 {range}

Cython: ----------------------------------------------------
 3 function calls in 2.352 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
 1 0.000 0.000 2.352 2.352 <string>:1(<module>)
 1 2.352 2.352 2.352 2.352 {Cythonteste.juntar_texto}
 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}

O Cython possui ainda outros recursos que este simples exemplo não aborda, vale a pena dar uma no site.

Anúncios

Mapeando Arvores Isoladas (ou qualquer outra feição pontual)

Em um trabalho recente o objetivo era fazer o levantamento de centenas de arvores em um determinado local. Por fim, além de um mapa, era necessário gerar um relatório contendo informações e foto de cada arvore.

Fazer todo o processo manualmente e elaborar o relatório em um documento comum do Word seria uma tarefa extremamente demorada e caso fosse necessário fazer alguma mudança, como no layout do documento por exemplo, o trabalho de alteração seria enorme. Assim resolvi automatizar o processo no seguinte fluxo de trabalho:

Coleta de dados -> Importar dados (Python) -> Banco de dados (PostgreSQL) -> Gerar Relatório (iReports) -> Gerar Mapa (Qgis)

O foco deste post é demonstrar a modelagem do banco de banco de dados e a importação dos dos dados, preparando o terreno para o uso do iReports e do Qgis, os quais não vou abordar aqui.

Este modelo foi utilizado para o caso das árvores, mas com pequenas alterações pode ser extrapolado para qualquer outro tipo de objeto pontual.

1) A primeira etapa (após a coleta de dados) é definir o modelo de banco de dados, após alguns rascunhos no papel cheguei ao modelo que segue:

Basicamente, cada árvore mapeada possui:
-Uma Espécie.
-Diversos Atributos.
-Apenas um ponto.
-Diversas fotos.

2)Codificar o modelo no SqlAlchemy.

Importar dados em um banco de dados relacional pode ser um processo complicado devido a presença de chaves externas, para facilitar a tarefa utilizei o ORM do SqlAlchemy e eliminei dois processos em um só:

1)O SqlAlchemy criou a estrutura do banco de dados para mim.

2)O banco de dados foi mapeado em objetos Python para que os dados sejam importados facilmente.

O código a seguir é a definição do modelo da imagem acima para ser interpretado pelo SqlAlchemy:

from sqlalchemy import create_engine
from sqlalchemy.ext.declarative import declarative_base
from sqlalchemy.orm import relationship, sessionmaker
from sqlalchemy import Column, Integer, String,Sequence,Boolean,ForeignKey,Numeric
from geoalchemy import *

#configurar o banco de dados aqui
engine = create_engine("postgresql+psycopg2://postgres:omoomo@localhost/postgis", echo=False)
Session=sessionmaker(bind=engine)
Base=declarative_base()

class Arvore(Base):
    __tablename__='arvores'
    __table_args__ = {'schema':'public'}
    id=Column(Integer,Sequence('arvores_id_seq'), primary_key=True)
    plaqueta=Column(Integer)
    especie_id = Column(Integer, ForeignKey('public.especies.id'))
    especie = relationship("Especie")
    foto = relationship("Foto")
    fitossanidade=Column(String)
    hfuste=Column(Numeric)
    dap=Column(Numeric)
    volume=Column(Numeric)
    ponto_id=Column(Integer, ForeignKey('public.pontos.id'))
    ponto=relationship("Ponto")

class Ponto(Base):
    __tablename__='pontos'
    __table_args__ = {'schema':'public'}
    id=Column(Integer, primary_key=True)
    the_geom = GeometryColumn(Point(dimension=2,srid=32723))

class Especie(Base):
    __tablename__='especies'
    __table_args__ = {'schema':'public'}
    id=Column(Integer,Sequence('especies_id_seq'), primary_key=True)
    especie=Column(String)
    nomepopular=Column(String)
    nativa=Column(Boolean)
    ameacada=Column(String)
    repouso=Column(String)

class Foto(Base):
    __tablename__='fotos'
    __table_args__ = {'schema':'public'}
    id=Column(Integer,Sequence('fotos_id_seq'), primary_key=True)
    arvore_id = Column(Integer, ForeignKey('public.arvores.id'))
    foto=Column(String)

GeometryDDL(Ponto.__table__)

if __name__=="__main__":
    Base.metadata.create_all(engine)

Não é necessário criar as tabelas no banco de dados, o comando

Base.metadata.create_all(engine)

cuida disto automaticament, e como este módulo será chamado diversas vezes mais para frente, este comando foi colocado dentro da condicional if __name__==”__main__” para evitar que quando for chamado ele tente criar as tabelas repetidamente.

3) A próxima etapa é importar os dados.

Para importar o arquivo gpx com os pontos utilizei a API do OGR (biblioteca que lê e escreve inúmeros formatos de arquivos de geoprocessamento).
Para importar os demais dados exportei as tabelas em formato de texto delimitado por tabulações e no Python lí os arquivos com a ajuda da biblioteca CSV.

#encoding:utf-8
#importa arquivos delimitados por tabulações

import model
import csv
from osgeo import ogr
from geoalchemy import WKTSpatialElement

session=model.Session()

#Arquivo de pontos GPX
pontos = ogr.Open( "./dados/pontos.gpx" )
camada=pontos.GetLayerByName('waypoints')
camada.ResetReading()
for ponto in camada:
    nome=ponto.GetField('name')
    geometria=ponto.GetGeometryRef()
    print geometria.ExportToWkt()
    try:
        nome=int(nome)
        itemponto=model.Ponto(id=nome,the_geom=WKTSpatialElement(geometria.ExportToWkt()))
        session.add(itemponto)
    except: pass

#arquivo de espécies
especies = csv.reader(open('./dados/especies.txt', 'rb'),dialect='excel-tab')
for especie in especies:
    itemespecie=model.Especie(id=float(especie[0].strip()),especie=especie[1])
    session.add(itemespecie)

#arquivo de arvores
arvores = csv.reader(open('./dados/arvorestab.txt', 'rb'),dialect='excel-tab')
listaplaquetas=[]
for arvore in arvores:
    try:hfuste=float(arvore[3].replace(',','.'))
    except:hfuste=0
    try:dap=float(arvore[4].replace(',','.'))
    except:dap=0
    if arvore[1]=='':especie_id = None
    else:especie_id=arvore[1]
    if arvore[5]=='':ponto_id=None
    else:ponto_id=arvore[5]
    itemarvore=model.Arvore(id=arvore[0],plaqueta=arvore[0],
                            especie_id=especie_id,fitossanidade=arvore[2],
                            hfuste=hfuste,dap=dap,ponto_id=ponto_id)
    listaplaquetas.append(arvore[0])
    session.add(itemarvore)

#fotos
for i in listaplaquetas:
    itemfoto=model.Foto(arvore_id=i,foto="\\fotos\\%s.jpg"%(i))
    session.add(itemfoto)

session.commit()

5) Gerar o Relatório.

Para gerar os campos no iReports é necessária uma SQL Query, eu uso o pgAdmin para facilitar o processo, crio parte da busca no editor visual e depois faço pequenas alterações no código. Na figura ao lado está a busca como fica no editor e abaixo o código gerado.

SELECT
  arvores.plaqueta,
  especies.especie,
  especies.nomepopular,
  especies.nativa,
  especies.ameacada,
  especies.repouso,
  arvores.fitossanidade,
  arvores.hfuste,
  arvores.dap,
  arvores.volume,
  ST_X(pontos.the_geom),
  ST_Y(pontos.the_geom),
  fotos.foto
FROM
  public.arvores,
  public.pontos,
  public.fotos,
  public.especies
WHERE
  arvores.especie_id = especies.id AND
  arvores.ponto_id = pontos.id AND
  fotos.arvore_id = arvores.id;

Para quem não conhece a ferramenta, o iReports é um gerador de relatórios com uma interface gráfica excelente, nele você adiciona os campos resultantes da busca SQL em um layout e com um comando ele gera o relatório em diversos formatos, incluindo PDF, Word e HTML.

De posse da busca SQL resta apenas desenhar o layout do relatório no iReports, mas este é um tema que vou abordar em um próximo post. Att.

Vetorização de Cartas Impressas no GRASS

O SIG GRASS possui uma excelente ferramenta para converter mapas e cartas impressas (e escaneadas) em linhas vetoriais.

O procedimento para a realização desta tarefa é o seguinte:

Primeiro é necessário filtrar o conteúdo a ser vetorizado, uma das formas de se fazer isto e abrindo o mapa em um editor de imagens como o GIMP, então use a ferramenta de selecionar por cor (select by colors) coloque alguma tolerância (15 deve ser suficiente) e selecione a cor preta. Inverta a seleção e apague todo o conteúdo que não esteja em preto. Salve como TIFF sem compressão.

Então no GRASS:

r.in.gdal – importar a imagem

r.thin – Afina as celulas não nulas que se tornarão linhas.
r.to.vect – Converte o raster em vetor.

Normalmente para um bom resultado é interessante fazer uma limpeza prévia da carta retirando legendas, sujeiras, etc.

Aqui vai um exemplo de uma carta já vetorizada:

Variações no Cálculo de Áreas

Algumas vezes que realizei o cálculo de áreas de polígonos e por algum motivo tive que verificar estas mesmas áreas em diferentes softwares me deparei com pequenas variações no resultado, alguns poucos metros quadrados de diferença para grandes áreas. Após me certificar que sistema de referência  utilizado era o mesmo em ambos os softwares resolvi pesquisar porque existe esta variação.

Obtive algumas respostas na lista ArcGis-Brasil e no site gis.stackexchange.com , basicamente existem diferentes fórmulas para calcular áreas, e diferentes softwares podem utilizar fórmulas diferentes e esta é a mais provável das causas desta variação.

Uma das formulas mais comuns é a de Simpson (http://www.spatialanalysisonline.com/).

Fast SQL Layer – Plugin para o QGIS

Para quem trabalha com banco de dados espaciais e camadas Postgis, utilizar os plugins “Add PostGIS/SpatiaLite Layer” e “RT Sql Layer” começa a ficar irritante pois toda a vez que quiser fazer uma mudança na camada, seja a menor que for, é necessário clicar em vários menus e digitar ou colar novamente o código.

É neste ponto que o plugin Fast SQL Layer otimiza o processo permitindo que camadas SQL sejam adicionadas rapidamente.

Continuar lendo