In [1]:
# Site selection analysis for new venture.
# organisation seeks housing for two functions in Randstad
# Requirements
''' 1. Function A: 1000-1500m2 as office/meeting room
2. Function B: 1000-3000m2 as logistics space (storage)
3. Maximum 30minutes travel time from Hoofddorp
4. Good logistical access - main traffic network
'''
Out[1]:
In [ ]:
In [2]:
# Import Libraries
import requests
import psycopg2, psycopg2.extras
import json
import math
import pandas as pd
import numpy as np
import shapely
from shapely import wkb
from shapely.geometry import box
import geopandas as gpd
import ipyleaflet
from ipyleaflet import *
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('bmh')
plt.rcParams['figure.figsize'] = (10.0, 10.0)
In [3]:
# Common Functions
def chunker(l, n):
"""Yield successive n-sized chunks from l."""
for i in range(0, len(l), n):
yield l[i:i + n]
def chunk_it(a, n):
b = chunker(a, n)
b = list(b)
return b
In [4]:
# GIS database connection
global gisdb
def postgres_test():
global gisdb
try:
cur = gisdb.cursor()
cur.execute('SELECT 1')
cur.close()
return gisdb
except:
gisdb = psycopg2.connect(host='',database="gis", user='p', password='p', cursor_factory=psycopg2.extras.RealDictCursor)
return gisdb
def showalltables():
cur = postgres_test().cursor()
cur.execute('SELECT * FROM pg_catalog.pg_tables')
data = cur.fetchall()
#data = showalltables()
tables = []
for i in data:
if i['schemaname'] == 'public':
tables.append(i['tablename'])
#print(i['tablename'])
cur.close()
return tables
def querydb(q):
cur = postgres_test().cursor()
cur.execute(q)
data = cur.fetchall()
cur.close()
return data
# Use rollback_db if you get any transaction related error
def rollback_db():
curs = postgres_test().cursor()
curs.execute("ROLLBACK")
cur.close()
gisdb.commit()
In [5]:
# Lists All tables in Database
showalltables()
Out[5]:
In [6]:
# Lists ALl gemmente
query = 'Select gm_code, gm_naam from gem_2017 ORDER BY gm_naam ASC LIMIT 10000'
querydb(query)
Out[6]:
In [ ]:
In [ ]:
In [7]:
# Variables
Hoofddorp_CentreCoordinates = [52.3061, 4.6907]
In [ ]:
In [8]:
# The number of monuments in Voorschoten
query = "Select count(*) as National_monuments from rijksmonumenten, (Select wkb_geometry as gm_boundary from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Haarlemmermeer') a WHERE ST_Intersects( wkb_geometry, a.gm_boundary)"
querydb(query)
Out[8]:
In [9]:
# The number of monuments in Haarlemmermeer and by type?
#query = "Select rijksmonumenten.* from rijksmonumenten, (Select wkb_geometry as gm_boundary from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Voorschoten') a WHERE ST_Intersects( wkb_geometry, a.gm_boundary) LIMIT 1000"
query = '''SELECT r.*, b.shape_area as builtarea
FROM ( Select rijksmonumenten.*
from rijksmonumenten,
(Select wkb_geometry as gm_boundary from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Haarlemmermeer') a
WHERE ST_Intersects( wkb_geometry, a.gm_boundary) LIMIT 1000
) r
LEFT JOIN kadastralekaartv3_bebouwing_poly b
ON ST_Within(r.wkb_geometry, b.wkb_geometry)
LIMIT 1000
'''
con = postgres_test()
df = gpd.GeoDataFrame.from_postgis(query, con, geom_col='wkb_geometry' )
print("Found Total", len(df), "Records")
df.groupby(['cbscategor']).size().reset_index(name='counts')
Out[9]:
In [10]:
df.plot(column='cbscategor', legend=True)
Out[10]:
In [11]:
# The number of monuments around Haarlemmermeer ? in meters 50000 = 50 km
#query = "Select rijksmonumenten.* from rijksmonumenten, (Select wkb_geometry as gm_boundary from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Voorschoten') a WHERE ST_Intersects( wkb_geometry, a.gm_boundary) LIMIT 1000"
query = '''SELECT r.*, b.shape_area as builtarea
FROM ( Select rijksmonumenten.*
from rijksmonumenten,
(Select wkb_geometry as gm_boundary from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Haarlemmermeer') a
WHERE ST_DWithin( wkb_geometry, a.gm_boundary, 5000) LIMIT 5000
) r
LEFT JOIN kadastralekaartv3_bebouwing_poly b
ON ST_Within(r.wkb_geometry, b.wkb_geometry)
'''
con = postgres_test()
df = gpd.GeoDataFrame.from_postgis(query, con, geom_col='wkb_geometry' )
print("Found Total", len(df), "Records")
df.groupby(['cbscategor']).size().reset_index(name='counts')
Out[11]:
In [12]:
#df.plot(column='cbscategor', legend=True)
df.plot(legend=True)
Out[12]:
In [13]:
# Get monuments based on travel time from point
def get_traveltime(startpoint, endpoints):
chunks = chunk_it(endpoints, 200)
durations = []
for chunk in chunks:
points = ";".join(chunk)
url = "https://live.renovatio.app/osrmapi/table/v1/car/" + startpoint + ";" + points + "?sources=0"
r = requests.get(url)
durations += r.json()['durations'][0][1:]
return durations
'''wgsdf = df.to_crs({'init': 'epsg:4326'})
points = list(wgsdf.wkb_geometry.x.astype(str) + ',' + wgsdf.wkb_geometry.y.astype(str))
traveltimes = get_traveltime('4.6907,52.3061', points)
wgsdf['traveltime'] = pd.Series(traveltimes)
wgsdf.plot(column='traveltime', scheme='fisher_jenks', legend=True)'''
Out[13]:
In [14]:
'''wgsdf'''
Out[14]:
In [15]:
# Filter Based on Travel time
'''drivetime = 30 #minutes
wgsdf = wgsdf[wgsdf['traveltime'] < drivetime*60]
wgsdf'''
Out[15]:
In [16]:
# Filter Based on AREA
'''builtarea_filter = (1000,2000)
wgsdf = wgsdf[(wgsdf['builtarea'] <= builtarea_filter[1]) & (wgsdf['builtarea'] >= builtarea_filter[0])]
wgsdf'''
Out[16]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [17]:
# Fishnet Algorithm
def fishnet(geometry, threshold):
bounds = geometry.total_bounds
xmin = int(bounds[0] // threshold)
xmax = int(bounds[2] // threshold)
ymin = int(bounds[1] // threshold)
ymax = int(bounds[3] // threshold)
ncols = int(xmax - xmin + 1)
nrows = int(ymax - ymin + 1)
result = []
for i in range(xmin, xmax+1):
for j in range(ymin, ymax+1):
b = box(i*threshold, j*threshold, (i+1)*threshold, (j+1)*threshold)
result.append(b)
return result
In [18]:
# Get Buffer From the Municipality
query = '''
Select wkb_geometry from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Haarlemmermeer'
'''
con = postgres_test()
mundf = gpd.GeoDataFrame.from_postgis(query, con, geom_col='wkb_geometry' )
mundf.plot()
Out[18]:
In [19]:
# Buffer
munBufferdf = mundf.buffer(20000)
munBufferdf.plot()
Out[19]:
In [20]:
# Create Fishnet
fishnet_grids = fishnet(munBufferdf, 1000)
crs = {'init': 'epsg:28992'}
fishnet_df = gpd.GeoDataFrame(fishnet_grids, crs=crs, geometry=fishnet_grids)
fishnet_df = fishnet_df.reset_index().drop(columns=[0])
fishnet_df['id'] = fishnet_df['index'] + 1000
print("Total", len(fishnet_df), "Grids found")
fishnet_df.plot(column='id')
Out[20]:
In [21]:
# Get Travel Time for Fishnets
fishnet_df_ = fishnet_df.to_crs({'init': 'epsg:4326'})
fishnet_df_['x'] = fishnet_df_['geometry'].centroid.x
fishnet_df_['y'] = fishnet_df_['geometry'].centroid.y
points = list(fishnet_df_.x.astype(str) + ',' + fishnet_df_.y.astype(str))
traveltimes = get_traveltime('4.6907,52.3061', points)
fishnet_df['traveltime'] = pd.Series(traveltimes)
fishnet_df.plot(column='traveltime', scheme='fisher_jenks', legend=True)
Out[21]:
In [ ]:
In [ ]:
In [22]:
# The number of monuments around Haarlemmermeer ? in meters 50000 = 50 km
#query = "Select rijksmonumenten.* from rijksmonumenten, (Select wkb_geometry as gm_boundary from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Voorschoten') a WHERE ST_Intersects( wkb_geometry, a.gm_boundary) LIMIT 1000"
query = '''SELECT r.*, b.shape_area as builtarea
FROM ( Select rijksmonumenten.*
from rijksmonumenten,
(Select wkb_geometry as gm_boundary from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Haarlemmermeer') a
WHERE ST_DWithin( wkb_geometry, a.gm_boundary, 20000) LIMIT 20000
) r
LEFT JOIN kadastralekaartv3_bebouwing_poly b
ON ST_Within(r.wkb_geometry, b.wkb_geometry)
'''
con = postgres_test()
df = gpd.GeoDataFrame.from_postgis(query, con, geom_col='wkb_geometry' )
print("Found Total", len(df), "Records")
df.groupby(['cbscategor']).size().reset_index(name='counts')
Out[22]:
In [23]:
df.plot()
Out[23]:
In [24]:
# Plot Both Layers
base = fishnet_df.to_crs({'init': 'epsg:28992'}).plot(column='traveltime', scheme='fisher_jenks', legend=True)
df.to_crs({'init': 'epsg:28992'}).plot(ax=base, marker='o', color='red', markersize=5, figsize=(12,12));
In [25]:
# filter grids based on travel time
drivetime = 30 #minutes
inDriveTime = fishnet_df[fishnet_df['traveltime'] < drivetime*60]
base = fishnet_df.to_crs({'init': 'epsg:28992'}).plot(column='traveltime', scheme='fisher_jenks', legend=True)
inDriveTime.to_crs({'init': 'epsg:28992'}).plot(ax=base, facecolor="none", edgecolor='red', lw=0.7)
Out[25]:
In [26]:
# Clip Monuments by inDriveTime
idtmon = df[df.intersects(inDriveTime.unary_union)]
print("Found", len(idtmon), "within", drivetime, "minutes drive from municipality")
idtmon.plot()
Out[26]:
In [27]:
# Plot All together
base = fishnet_df.to_crs({'init': 'epsg:28992'}).plot(column='traveltime', scheme='fisher_jenks', legend=True)
base = inDriveTime.to_crs({'init': 'epsg:28992'}).plot(ax=base, facecolor="none", edgecolor='red', lw=0.7)
idtmon.to_crs({'init': 'epsg:28992'}).plot(ax=base, marker='o', color='white', markersize=5);
In [28]:
list = idtmon[idtmon['cbscategor'] == 'Agrarische gebouwen']
list
Out[28]:
In [29]:
idtmon
Out[29]:
In [ ]:
In [ ]:
In [30]:
# Get Transport Data
query = '''
Select rijstroken.* from rijstroken,
(Select wkb_geometry as gm_boundary from bestuurlijkegrenzen_gemeentegrenzen WHERE gemeentenaam='Haarlemmermeer') a
WHERE ST_DWithin( wkb_geometry, a.gm_boundary, 30000) LIMIT 20000
'''
con = postgres_test()
tpdf = gpd.GeoDataFrame.from_postgis(query, con, geom_col='wkb_geometry' )
tpdf.plot() # Highway data right now
Out[30]:
In [31]:
# Plot Highways with Municpal Boundaries And buffer
tpbase = mundf.plot(facecolor="none", edgecolor='red', lw=0.7)
tpbase = munBufferdf.plot(ax=tpbase, facecolor="none", edgecolor='red', lw=0.7)
tpdf.plot(ax=tpbase)
Out[31]:
In [32]:
# Convert Highways to single geometry
tpdf['a'] = 0
tp_dis_df = tpdf.dissolve(by='a')
tpbase = mundf.plot(facecolor="none", edgecolor='red', lw=0.7)
tpbase = munBufferdf.plot(ax=tpbase, facecolor="none", edgecolor='red', lw=0.7)
tp_dis_df.plot(ax=tpbase)
Out[32]:
In [ ]:
In [33]:
# Buffer the Highways 100M
tpbase = mundf.plot(facecolor="none", edgecolor='black', lw=0.7)
tpbase = munBufferdf.plot(ax=tpbase, facecolor="none", edgecolor='red', lw=0.7)
tp500df = gpd.GeoDataFrame(gpd.GeoSeries(tpdf.buffer(100).unary_union))
tp500df.columns = ['geometry']
tpbase = tpdf.plot(ax=tpbase)
tpbase = tp500df.plot(ax=tpbase, facecolor="none", edgecolor='blue', lw=0.7)
idtmon500tp = idtmon[idtmon.intersects(tp500df.unary_union)]
print("Found", len(idtmon500tp), "within 100 Mters buffer from Highway")
idtmon500tp.plot(ax=tpbase, marker='o', color='red', markersize=5);
In [34]:
idtmon500tp.columns
Out[34]:
In [35]:
List100 = idtmon500tp[['rijksmonnr', 'cbscategor', 'eindbouwjr','builtarea','oorsprfunc','subcatoms', 'straat','huisnummer','toevoeging','postcode','gemeente']]
List100.to_excel('output.xlsx')
In [36]:
# Buffer the Highways 200M
tpbase = mundf.plot(facecolor="none", edgecolor='black', lw=0.7)
tpbase = munBufferdf.plot(ax=tpbase, facecolor="none", edgecolor='red', lw=0.7)
tp500df = gpd.GeoDataFrame(gpd.GeoSeries(tpdf.buffer(200).unary_union))
tp500df.columns = ['geometry']
tpbase = tpdf.plot(ax=tpbase)
tpbase = tp500df.plot(ax=tpbase, facecolor="none", edgecolor='blue', lw=0.7)
idtmon500tp = idtmon[idtmon.intersects(tp500df.unary_union)]
print("Found", len(idtmon500tp), "within 200 Mters buffer from Highway")
idtmon500tp.plot(ax=tpbase, marker='o', color='red', markersize=5);
In [37]:
# Buffer the Highways 500M
tpbase = mundf.plot(facecolor="none", edgecolor='black', lw=0.7)
tpbase = munBufferdf.plot(ax=tpbase, facecolor="none", edgecolor='red', lw=0.7)
tp500df = gpd.GeoDataFrame(gpd.GeoSeries(tpdf.buffer(500).unary_union))
tp500df.columns = ['geometry']
tpbase = tpdf.plot(ax=tpbase)
tpbase = tp500df.plot(ax=tpbase, facecolor="none", edgecolor='blue', lw=0.7)
idtmon500tp = idtmon[idtmon.intersects(tp500df.unary_union)]
print("Found", len(idtmon500tp), "within 500 Mters buffer from Highway")
idtmon500tp.plot(ax=tpbase, marker='o', color='red', markersize=5);
In [38]:
# Buffer the Highways 1KM
tpbase = mundf.plot(facecolor="none", edgecolor='black', lw=0.7)
tpbase = munBufferdf.plot(ax=tpbase, facecolor="none", edgecolor='red', lw=0.7)
tp500df = gpd.GeoDataFrame(gpd.GeoSeries(tpdf.buffer(1000).unary_union))
tp500df.columns = ['geometry']
tpbase = tpdf.plot(ax=tpbase)
tpbase = tp500df.plot(ax=tpbase, facecolor="none", edgecolor='blue', lw=0.7)
idtmon500tp = idtmon[idtmon.intersects(tp500df.unary_union)]
print("Found", len(idtmon500tp), "within 500 Mters buffer from Highway")
idtmon500tp.plot(ax=tpbase, marker='o', color='red', markersize=5);
In [39]:
# Buffer the Highways 2KM
tpbase = mundf.plot(facecolor="none", edgecolor='black', lw=0.7)
tpbase = munBufferdf.plot(ax=tpbase, facecolor="none", edgecolor='red', lw=0.7)
tp500df = gpd.GeoDataFrame(gpd.GeoSeries(tpdf.buffer(2000).unary_union))
tp500df.columns = ['geometry']
tpbase = tpdf.plot(ax=tpbase)
tpbase = tp500df.plot(ax=tpbase, facecolor="none", edgecolor='blue', lw=0.7)
idtmon500tp = idtmon[idtmon.intersects(tp500df.unary_union)]
print("Found", len(idtmon500tp), "within 500 Mters buffer from Highway")
idtmon500tp.plot(ax=tpbase, marker='o', color='red', markersize=5);
In [ ]:
In [ ]: