What spiega is about
- Knowledge sharing
- Portfolio
geo transformation and queries for GIS
For many years I worked with GIS software to run analytics on spatial data. Here is a collection of different operations I used across different projects which I wanted to share to people who might need to run similar analysis.
Geo operations
Materials:
*
Points in polygon
Library
We need the libraries:
import os, sys, gzip, random, csv, json, datetime, re
import numpy as np
import pandas as pd
import scipy as sp
import matplotlib.pyplot as plt
import geopandas as gpd
from scipy.spatial import cKDTree
from scipy import inf
import shapely as sh
baseDir = "~/raw/"
Data
poi = pd.read_csv(baseDir + "poi.csv")
#x,y,name,region,bla bla
Points in polygon
Coordinates in region:
region = gpd.GeoDataFrame.from_file(baseDir + "gis/geo/bundesland.shp")
region.index = region['GEN']
region = region['geometry']
pL = poi[['x','y']].apply(lambda x: sh.geometry.Point(x[0],x[1]),axis=1)
pnts = gpd.GeoDataFrame(geometry=pL)
pnts = pnts.assign(**{key: pnts.within(geom) for key, geom in region.items()})
for i in pnts.columns[1:]:
poi.loc[pnts[i],"region"] = i
Clustering
Cluster coordinates within a distance
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.cluster.hierarchy import fcluster
Z = linkage(poi[['x','y']], 'ward')
gradMeter = 111122.19769899677
max_d = 1500./gradMeter
poi.loc[:,'id_zone'] = fcluster(Z,max_d,criterion='distance')
Angle
Calculate orientation of a segment, orthogonal distance from segment to reference point, calculates chirality
dist1 = np.sqrt((nodeL['x1'] - nodeL['x_poi'])**2 + (nodeL['y1'] - nodeL['y_poi'])**2)
dist2 = np.sqrt((nodeL['x2'] - nodeL['x_poi'])**2 + (nodeL['y2'] - nodeL['y_poi'])**2)
nodeL.loc[:,'dist'] = np.min([dist1,dist2],axis=0)
nodeL.loc[:,"orth_dist"] = np.abs((nodeL['x2']-nodeL['x1'])*(nodeL['y1']-nodeL['y_poi'])-(nodeL['y2']-nodeL['y1'])*(nodeL['x1']-nodeL['x_poi']))
nodeL.loc[:,"orth_dist"] = nodeL["orth_dist"]/(np.abs((nodeL['x2']-nodeL['x1'])) + np.abs((nodeL['y2']-nodeL['y1'])))
nodeL.loc[:,"dist"] = nodeL['orth_dist']*nodeL['dist']
v1 = [nodeL['x1'] - metr['deCenter'][0],nodeL['y1'] - metr['deCenter'][1]]
v2 = [nodeL['x2'] - nodeL['x1'],nodeL['y2'] - nodeL['y1']]
crossP = v1[0]*v2[1] - v2[0]*v1[1]
nodeL.loc[:,'chirality_v'] = 1.*(crossP > 0.)
Chirality
Calculates the chirality between two angles
poi.loc[:,'angle'] = np.arctan2((poi['y']-poi['y_mot']),(poi['x']-poi['x_mot']))*180./np.pi
poi.loc[:,'angle'] = - poi.loc[:,'angle']
poi.loc[:,'tang'] = np.arctan2(metr['deCenter'][1]-poi['y'],metr['deCenter'][0]-poi['x'])*180./np.pi
poi.loc[:,'tang'] = 90. - poi.loc[:,'tang']
poi.loc[poi['tang']>180.,'tang'] -= 180.
t = np.abs(poi['tang']-poi['angle'])
t[t>180.] = 360.-t
poi.loc[:,'chirality'] = 1*(t>90)
Tangent point
Calculates the tangent point on a line from a reference point
from shapely.ops import split, snap
from shapely import geometry, ops
motG = gpd.GeoDataFrame.from_file(baseDir + "gis/geo/motorway.shp")
motG = motG[motG['geometry'].apply(lambda x: x.is_valid).values]
line = motG.geometry.unary_union
for i,poii in poi.iterrows():
p = geometry.Point(poi.loc[i][['x','y']])
neip = line.interpolate(line.project(p))
#snap(coords, line, tolerance)
poi.loc[i,"x_mot"] = neip.x
poi.loc[i,"y_mot"] = neip.y
Contains
Which polygon contains a point
import shapely.speedups
shapely.speedups.enable()
densG = gpd.GeoDataFrame.from_file(baseDir + "gis/geo/pop_dens_2km.shp")
g = densG['geometry'][0]
p = [sha.geometry.Point(x,y) for x,y in zip(poi['x'],poi['y'])]
poiG = gpd.GeoDataFrame(p,columns=["geometry"])
Polygon to edges
Dissolve a polygon into its edges
tileG = gpd.read_file(baseDir + "gis/tank/tileList.geojson")
tileG.loc[:,'sum'] = tileG.loc[:,'north_in'] + tileG.loc[:,'south_in'] + tileG.loc[:,'east_in'] + tileG.loc[:,'west_in']
dirg = tileG[['tile_id','sum','north_in','south_in','east_in','west_in','north_out','south_out','east_out','west_out']].groupby(['tile_id']).agg(sum)
dirg = dirg.reset_index()
tileL = tileG[['tile_id','col_id','row_id','geometry']].groupby(['tile_id']).head(1)
dirg = pd.merge(dirg,tileL,left_on="tile_id",right_on="tile_id",how="left")
dirg = gpd.GeoDataFrame(dirg)
with open(baseDir + "gis/tank/junction_tile.geojson","w") as fo:
fo.write(dirg.to_json())
dirl = gpd.GeoDataFrame(columns=["in","out","dir","geometry"])
for i,a in dirg.iterrows():
l = a['geometry'].boundary
ll = LineString([(l.xy[0][0],l.xy[1][0]),(l.xy[0][1],l.xy[1][1])])
dirl.loc[str(a['tile_id']) + 'a'] = [a['east_in'],a['east_out'],"e",ll]
ll = LineString([(l.xy[0][1],l.xy[1][1]),(l.xy[0][2],l.xy[1][2])])
dirl.loc[str(a['tile_id']) + 'b'] = [a['north_in'],a['north_out'],"n",ll]
ll = LineString([(l.xy[0][2],l.xy[1][2]),(l.xy[0][3],l.xy[1][3])])
dirl.loc[str(a['tile_id']) + 'c'] = [a['west_in'],a['west_out'],"w",ll]
ll = LineString([(l.xy[0][3],l.xy[1][3]),(l.xy[0][4],l.xy[1][4])])
dirl.loc[str(a['tile_id']) + 'd'] = [a['south_in'],a['south_out'],"s",ll]
dirl = gpd.GeoDataFrame(dirl)
with open(baseDir + "gis/tank/junction_edge.geojson","w") as fo:
fo.write(dirl.to_json())
Spectral clustering
from sklearn.cluster import SpectralClustering
from sklearn.cluster import KMeans
mat = np.matrix([[1.,.1,.6,.4],[.1,1.,.1,.2],[.6,.1,1.,.7],[.4,.2,.7,1.]])
print(SpectralClustering(2).fit_predict(mat))
eigen_values, eigen_vectors = np.linalg.eigh(mat)
print(KMeans(n_clusters=2, init='k-means++').fit_predict(eigen_vectors[:, 2:4]))
from sklearn.cluster import DBSCAN
DBSCAN(min_samples=1).fit_predict(mat)
Intersect lines and areas
junct = gpd.read_file(baseDir + "junction_area.geojson")
dirc = gpd.read_file(baseDir + "count_dir.geojson")
dirA = [np.arctan(0.),np.arctan(np.pi/2.),np.arctan(np.pi),np.arctan(3.*np.pi/2.),np.arctan(2.*np.pi)]
def getAng(dx,dy):
ang = int(np.arctan2(dy,dx)*2./np.pi + 0.5)
cordD = ["east","north","west","south"]
return cordD[ang], cordD[abs(2-ang)]
def getDir(dtx,dty):
cordD = [("east","west"),("north","south"),("west","east"),("south","north")]
ang = 1
if(dtx < 0):
ang = 0
elif(dtx > 0):
ang = 2
elif(dty < 0):
ang = 3
return cordD[ang]
cordD = [("east","west"),("north","south"),("west","east"),("south","north")]
exits = gpd.read_file(baseDir + "/motorway_exit_axes.geojson")
fluxC = pd.DataFrame(index=range(0,4*24))
exits.loc[:,'exit'] = 0
exits.loc[:,'enter'] = 0
for i,ex in exits.iterrows():
l = ex['geometry']
inTile = [a.contains(Point(l.xy[0][0],l.xy[1][0])) for a in dirc['geometry']]
outTile = [a.contains(Point(l.xy[0][1],l.xy[1][1])) for a in dirc['geometry']]
dircI = dirc.loc[inTile]
dircO = dirc.loc[outTile]
# ang1, ang2 = getAng(l.xy[1][1] - l.xy[1][0],l.xy[0][1] - l.xy[0][0])
dtx = dircI.iloc[0]['col_id'] - dircO.iloc[0]['col_id']
dty = dircI.iloc[0]['row_id'] - dircO.iloc[0]['row_id']
ang1, ang2 = getDir(dtx,dty)
Resampling
densG = gpd.GeoDataFrame.from_file(baseDir + "gis/geo/pop_density.shp")
centL = densG['geometry'].apply(lambda x: x.centroid)
densG.loc[:,"hash"] = centL.apply(lambda x: geohash.encode(x.xy[0][0],x.xy[1][0],precision=5))
def clampF(x):
return pd.Series({"pop_dens":x['Einwohner'].sum()
,"flat_dens":x['Wohnfl_Bew'].sum()
,"foreign":x['Auslaender'].sum()
,"women":x['Frauen_A'].sum()
,"young":x['unter18_A'].sum()
,"geometry":cascaded_union(x['geometry'])
,"household":x['HHGroesse_'].sum()
,"n":len(x['Flaeche'])
})
densG = densG.groupby("hash").apply(clampF).reset_index()
densG.loc[:,'geometry'] = densG['geometry'].apply(lambda f: f.convex_hull)
for i in ['pop_dens','flat_dens','foreign','women','young','household']:
densG.loc[:,i] = densG[i]/densG['n']
densG = gpd.GeoDataFrame(densG)
densG.to_file(baseDir + "gis/geo/pop_dens_2km.shp")
Mongo db
Library
import os, sys, gzip, random, csv, json, datetime, re
import numpy as np
import pandas as pd
import scipy as sp
import matplotlib.pyplot as plt
import geopandas as gpd
from scipy.spatial import cKDTree
from scipy import inf
import shapely as sh
import pymongo
baseDir = "~/raw/"
We initiate the client
with open(baseDir + '/credenza/geomadi.json') as f:
cred = json.load(f)
with open(baseDir + '/raw/metrics.json') as f:
metr = json.load(f)['metrics']
client = pymongo.MongoClient(cred['mongo']['address'],cred['mongo']['port'])
coll = client["index_name"]["collection_name"]
Near sphere
Returns all points within a distance
neiDist = 200.
nodeL = []
for i,poii in poi.iterrows():
poii = poi.loc[i]
poi_coord = [x for x in poii.ix[['x','y']]]
neiN = coll.find({'loc':{'$nearSphere':{'$geometry':{'type':"Point",'coordinates':poi_coord},'$minDistance':0,'$maxDistance':neiDist}}})
nodeId = []
for neii in neiN:
nodeL.append({'id_poi':poii['id_poi'],'src':neii['src'],'trg':neii['trg'],"maxspeed":neii['maxspeed'],'street':neii['highway']
,"x_poi":poii['x'],"y_poi":poii['y']
})
Intersects
Take all locations inside polygons
motG = gpd.GeoDataFrame.from_file(baseDir + "gis/geo/motorway_area.shp")
cellL = []
for g in np.array(motG['geometry'][0]):
c = g.exterior.coords.xy
c1 = [[x,y] for x,y in zip(c[0],c[1])]
neiN = coll.find({'geom':{'$geoIntersects':{'$geometry':{'type':"Polygon",'coordinates':[c1]}}}})
neii = neiN[0]
for neii in neiN:
cellL.append({"cilac":str(neii['cell_ci']) + '-' + str(neii['cell_lac'])})
cellL = pd.DataFrame(cellL)
Filtering
Filtering by list
coll = client["tdg_infra"]["infrastructure"]
poi = pd.read_csv(baseDir + "raw/tr_cilac_sel1.csv")
colL = list(poi.columns)
colL[0] = 'domcell'
poi.columns = colL
poi.loc[:,'ci'] = [re.sub("-.*","",x) for x in poi['domcell']]
poi.loc[:,'lac'] = [re.sub(".*-","",x) for x in poi['domcell']]
queryL = []
for i,p in poi.iterrows():
queryL.append({"cell_ci":p['ci']})
queryL.append({"cell_lac":p['lac']})
Boundary Box
neiN = coll.find({'loc':{'$geoWithin':{'$box':[ [BBox[0],BBox[2]],[BBox[1],BBox[3]] ]}}})
neo4j
from neo4j.v1 import GraphDatabase, basic_auth
driver = GraphDatabase.driver("bolt://localhost:7687", auth=basic_auth("neo4j", "neo4j"))
session = driver.session()
session.run("CREATE (a:Person {name: {name}, title: {title}})",
{"name": "Arthur", "title": "King"})
result = session.run("MATCH (a:Person) WHERE a.name = {name} "
"RETURN a.name AS name, a.title AS title",
{"name": "Arthur"})
for record in result:
print("%s %s" % (record["title"], record["name"]))
session.close()
from py2neo import Graph, Path
graph = Graph()
tx = graph.cypher.begin()
for name in ["Alice", "Bob", "Carol"]:
tx.append("CREATE (person:Person {name:{name}}) RETURN person", name=name)
alice, bob, carol = [result.one for result in tx.commit()]
friends = Path(alice, "KNOWS", bob, "KNOWS", carol)
graph.create(friends)
from neomodel import StructuredNode, StringProperty, RelationshipTo, RelationshipFrom, config
config.DATABASE_URL = 'bolt://neo4j:test@localhost:7687'
class Book(StructuredNode):
title = StringProperty(unique_index=True)
author = RelationshipTo('Author', 'AUTHOR')
class Author(StructuredNode):
name = StringProperty(unique_index=True)
books = RelationshipFrom('Book', 'AUTHOR')
harry_potter = Book(title='Harry potter and the..').save()
rowling = Author(name='J. K. Rowling').save()
harry_potter.author.connect(rowling)
Network
Library
import osmnx as ox
import networkx as nx
Graph:
graph1 = ox.load_graphml(filename="germany_split_motorway_motorwaylink.graphml")
graph2 = ox.load_graphml( filename="germany_split_trunk_trunk_link.graphml")
graph3 = ox.load_graphml( filename="germany_split_primlink.graphml")
graph4 = ox.load_graphml( filename="germany_split_prim.graphml")
graph5 = ox.load_graphml( filename="germany_split_seclink.graphml")
Compose:
c_graphs = [graph1, graph2, graph3, graph4, graph5]
composed_G = nx.compose_all(c_graphs)
Simplify:
simp_G = ox.simplify_graph(composed_G)
connected_G = max(nx.strongly_connected_component_subgraphs(simp_G), key=len)
graph_proj = ox.project_graph(connected_G)
ox.save_graphml(graph_proj, filename="allGermany_allstreetsUntilSec_proj.graphml")
Create network
import networkx as nx
G=nx.Graph()
G.add_node(1)
G.add_nodes_from([2,3])
H=nx.path_graph(10)
G.add_nodes_from(H)
G.add_node(H)
G.add_edge(1,2)
e=(2,3)
G.add_edge(*e)
G.add_edges_from([(1,2),(1,3)])
G.add_edges_from(H.edges())
G.remove_node(H)
G.clear()
G.add_edges_from([(1,2),(1,3)])
G.add_node(1)
G.add_edge(1,2)
G.add_node("spam")
G.add_nodes_from("spam")
G.number_o../../f/f_nodes()
G.number_o../../f/f_edges()
G.nodes()
G.edges()
G.neighbors(1)
G.remove_nodes_from("spam")
G.nodes()
G.remove_edge(1,3)
H=nx.DiGraph(G)
H.edges()
edgelist=[(0,1),(1,2),(2,3)]
H=nx.Graph(edgelist)
G.add_edge(1,3)
G[1][3]['color']='blue'
FG=nx.Graph()
FG.add_weighted_edges_from([(1,2,0.125),(1,3,0.75),(2,4,1.2),(3,4,0.375)])
for n,nbrs in FG.adjacency_iter():
for nbr,eattr in nbrs.items():
data=eattr['weight']
if data<0.5: print('(%d, %d, %.3f)' % (n,nbr,data))
for (u,v,d) in FG.edges(data='weight'):
if d<0.5: print('(%d, %d, %.3f)'%(n,nbr,d))
G = nx.Graph(day="Friday")
G.graph
G.graph['day']='Monday'
G.graph
import matplotlib.pyplot as plt
nx.draw(G)
nx.draw_random(G)
nx.draw_circular(G)
nx.draw_spectral(G)
plt.show()