Mineração de Dados Geoespaciais - Anexo I¶
Fernando Ferraz Ribeiro
Beth Leite Soares
O objetivo deste notebook é juntar dados de um shapefile com um mapa do Open Street Maps(OSM)
Testando o ambiente de trabalho (conda environment)¶
Mostrando o caminho que o servidor Jupyter carrega o interpretador Python (python.exe) e carregas o anbiente (\envs)
import os
import sys
print(sys.executable)
pathFix = sys.prefix
print(pathFix)
Lidando com erros do ambiente (WINDOWS ONLY)¶
Algumas bibliotecas importadas apresentam um erro apenas na plataforma Windows.
Caso esteja rodando no windows e tenha problemas com o sistema de coordenadas de referência (crs), rode a linha de comando abaixo:
pathFix = pathFix .replace('\\', '/')
print(pathFix)
os.environ["PROJ_LIB"] = pathFix + "/Library/share"
mais informações no Post sobre o Erro
Importando Bibliotecas¶
# Biblioteca basica de programação cientÃfica em python
import numpy as np
# biblioteca de análise de dados
import pandas as pd
# biblioteca de gráficos
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors
# Bibliotecas geopandas
import geopandas as gpd
# biblioteca de redes complexas
import networkx as nx
# biblioteca para acessar dados do Open sreet maps
import osmnx as ox
# ler dados remotamente da internet
%matplotlib inline
ox.config(use_cache=True, log_console=True)
ox.__version__
Criando coordenadas de recorte¶
# coordenada inicial x
xSC = 555000
# variação da coordenada x
deltaX = 2000
# coordenada inicial y
ySC = 8570000
# variação da coordenada y
deltaY = 2000
Limites da importaçâo¶
from shapely.geometry import Polygon
recorte = gpd.GeoSeries([ Polygon([(xSC,ySC), (xSC + deltaX , ySC), (xSC + deltaX, ySC + deltaY ), (xSC, ySC + deltaY )]) ])
A geometria criada pela linha de comando acima, embora tenha as corrdenadas relativas ao sistema de projeção Sigras 2000, não tem nenhuma informação georreferenciada. É preciso informar qual o sistema de coordinadas de referência utilizado ( coordenates reference sistem - crs). O Bloco de código abaxo informa que as coordenadas do recorte devem ser tratadas com o sistema Sigras 2000. com unidades em metro.
# colocando em coordenaadas SIRGAS 2000
recorte.crs = {'proj': 'utm', 'zone': 24, 'south': True, 'ellps': 'aust_SA', 'units': 'm', 'no_defs': True}
Carregando shapes¶
# importando Shape dos bairros - polÃgonos
bairros = gpd.read_file('../shapefiles/BaseSSA/Limites/bairros_fim.shp', bbox = recorte )
# importando Shape das edificaçõs - polilinhas
edf = gpd.read_file('../shapefiles/BaseSSA/edificacoes_polyline.shp', bbox = recorte)
# importando shape de pontos
edf_pt = gpd.read_file('../shapefiles/BaseSSA/edificacoes_point.shp', bbox = recorte)
Mostrando os dados dos shapes¶
dados do shape bairros¶
bairros.head()
Dados do shape edificacoes¶
edf.head()
edf_pt.head()
Plotando imagens¶
limites do recorte¶
fig, layers = plt.subplots(figsize=(10,10)
#,dpi=30
)
# pltando a imagem de fundo
recorte.plot(ax = layers, color = 'blue', alpha = 1)
# Plotando os bairros
# layers - quadro o qual poderá ser alterado os valores de xmin, xmax, tÃtulo e etc
# Plotando os bairros
bairros.plot(ax= layers, color='yellowgreen', edgecolor='black', linestyle="--", lw= 1.5, alpha=1)
# Plotando as edificações com a cor vermelha
edf.plot(ax= layers, color='darkred', lw= 0.7)
try:
edf_pt.plot(ax= layers, color='orange', lw= 0.7, )
except:
pass
# titulo da figura
fig.suptitle('Imagem Teste', fontsize=16)
# limites do gráfico
layers.set_xlim(xSC, xSC + deltaX )
layers.set_ylim(ySC, ySC + deltaY )
#fig.show()
Lendo Arquivos do banco de dados espaciais online Open Street maps¶
Através de arquivo baixado manualmente pelo site¶
gFile = ox.graph_from_file('../OSM/map.osm')
gFile
fig, ax = ox.plot_graph(gFile)
plt.tight_layout()
Baixando arquivos por coordenadas limite¶
O open street maps trabalha com coordenadas WGS84, definidas pelo código epsg 4326. O sistema WGS84 utiliza coordenadas em graus de latitude e longitude.
recorte.bounds
Para utilizar a mesma geometria limite da utilizada para a importação do Shapefile, é preciso converter o sistema de coordenadas de referência de Sigras 2000 para WGS84.
recorte_LL = recorte.to_crs({'init': 'epsg:4326'})
recorte_LL.bounds
O elemento de Ãndice 0 da coluna geometry é um polÃgono da biblioteca shapely.
type(recorte_LL.geometry[0])
segundo a documentação do comando ox.graph_from_polygon ele recebe como priméiro parâmetro um polÃgono ou multi-polÃgono da biblioteca citada.
O download do multi-grafo depende da conexão com a internet
gLatLon = ox.graph_from_polygon(
recorte_LL.geometry[0]
, network_type = 'all_private'
, truncate_by_edge = True
, retain_all = True
)
gLatLon
ox.plot_graph(ox.project_graph(gLatLon))
Através do endereço (geocoding)¶
gCode = ox.graph_from_place('Salvador, Bahia, Brazil', network_type="drive")
O download do multi-grafo depende da conexão com a internet
gCode
ox.plot_graph(ox.project_graph(gCode))
Centralidade é uma medida da importância de um vértice em um grafo. Os critérios de medida da importância variam. Os quatro mais utilizados são:
1. grau
2. intermediação
3. proximidade
4. vetor próprio
O código abaixo calcula a centralidade de proximidade da malha viária da localidade gravada no grafo do OSM
Antenção!¶
O custo computacional do cálculo da centralidade cresce exponencialmente de acordo com o número de vértices do grafo(a execução do código abaixo pode demorar muitas horas)
# edge closeness centrality: convert graph to line graph so edges become nodes and vice versa
edge_centrality = nx.closeness_centrality(nx.line_graph(gCode))
O código abaixo desenha um gráfico mostrando, em escala de cores a varição da medida de centralidade da malha viária.
# list of edge values for the orginal graph
ev = [edge_centrality[edge + (0,)] for edge in gCode.edges()]
# color scale converted to list of colors for graph edges
norm = colors.Normalize(vmin=min(ev)*0.8, vmax=max(ev))
cmap = cm.ScalarMappable(norm=norm, cmap=cm.inferno)
ec = [cmap.to_rgba(cl) for cl in ev]
# color the edges in the original graph with closeness centralities in the line graph
fig, ax = ox.plot_graph(gCode, bgcolor=(0.1,0.3,0.0,0.0), axis_off=True, node_size=0,
edge_color=ec, edge_linewidth=1.5, edge_alpha=1)
Juntando o arquivo do OSM com o Arquivo SHP¶
Trasnformando de Grafo OSM para geopandas Data frame¶
osm_pontos, osm_linhas = ox.save_load.graph_to_gdfs(gLatLon)
osm_linhas.head()
# Mudando o sistema de coordenadas de referência
# Para o Shape de Pontos
osm_pontos.to_crs(bairros.crs, inplace = True)
# Para o Shape de Bairros
osm_linhas.to_crs(bairros.crs, inplace = True)
Odenando as colunas do shape das linhas
# Ordenando as colunas das linhas para que a coluna de geometria fique por último
osm_linhas_cols = osm_linhas.columns.tolist()
osm_linhas_cols = [x for x in osm_linhas_cols if x != 'geometry']
osm_linhas_cols.append('geometry')
osm_linhas = osm_linhas[osm_linhas_cols]
osm_linhas.head()
osm_pontos.head()
Estatística descritivas do DF osm_linhas¶
osm_linhas.describe(include='all')
fig2, layers2 = plt.subplots(figsize=(10,10)
#,dpi=30
)
recorte.plot(ax = layers2, color = 'blue', alpha = 1)
# Plotando os bairros
# Limites preto
# alpha - transparência
# layers - quadro o qual poderá ser alterado os valores de xmin, xmax, tÃtulo e etc
bairros.plot(ax= layers2, color='yellowgreen', edgecolor='black', linestyle="--", lw= 1.5, alpha=1)
# Plotando as edificações com a cor vermelha
edf.plot(ax= layers2, color='darkred', lw= 0.7)
try:
edf_pt.plot(ax= layers2, color='orange', lw= 0.7, )
except:
pass
# importados do OSM
osm_linhas.plot(ax= layers2, color='teal', lw= 0.8, )
osm_pontos.plot(ax= layers2, color='indigo', lw= 0.4, )
# titulo da figura
fig2.suptitle('Imagem Teste', fontsize=16)
# limites do gráfico
layers2.set_xlim(xSC, xSC + deltaX )
layers2.set_ylim(ySC, ySC + deltaY )
osm_linhas.dtypes
Devido a incompatibilidades entre as estruturas de dados do shapefile e do OSM algumas colunas precisam ser tratadas
osm_linhas.oneway = osm_linhas.oneway.astype(int)
osm_linhas.key = osm_linhas.key.astype(int)
osm_linhas.name = osm_linhas.name.astype(str)
osm_linhas.bridge = osm_linhas.bridge.astype(str)
osm_select = osm_linhas.filter( items= [ 'name', 'oneway','key','length', 'bridge' ,'geometry' ])
osm_select.to_file('./output/teste01.shp', driver='ESRI Shapefile')