# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # San Jose Sidewalks Data Processing Documentation
#
# ## Summary
# Use Docker for POSTGIS, run some queries (I use Python Geopandas to get the data into POSTGIS and document the queries), validate and export using QGIS.
#
# There is probably a way easier way to do this is ESRI. I know.

# ### Get PostGIS (i'm using a Mac with docker).
#
# I'm using this instance:
# https://hub.docker.com/r/kartoza/postgis/
#
# Run Commands
# ```docker pull kartoza/postgis```
#
# ```docker run --name "postgis" -p 5432:5432 -d -t kartoza/postgis```
#
# username, password, both are docker
#
# #### Postgis Version
# POSTGIS="2.3.2 r15302" GEOS="3.4.2-CAPI-1.8.2 r3921" SFCGAL="1.2.2" 
# PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" 
# LIBXML="2.9.1" LIBJSON="0.11.99" TOPOLOGY RASTER

# +
import geopandas as gpd
import pandas as pd
import sqlalchemy as sal
from sqlalchemy import create_engine

# Connect to local docker POSTGIS instance
engine = sal.create_engine('postgresql://docker:docker@0.0.0.0/gis', echo=True)
conn = engine.connect()

# +
## Read in wgs 84 sidewalks file.
# I couldn't get geopandas to change the crs from the original data
# correctly so I used ESRI to do it. This is the output.
# There original data looks something like Basemap2_201707180931079098/SidewalkArea.shp'
sidewalks = gpd.read_file('/data/SJ_Sidewalks.shp')

## Load data into the POSTGIS server
# https://stackoverflow.com/questions/38361336/write-geodataframe-into-sql-database

# Function to generate WKB hex
def wkb_hexer(line):
    return line.wkb_hex

gdf = sidewalks
# Convert `'geom'` column in GeoDataFrame `gdf` to hex
# Note that following this step, the GeoDataFrame is just a regular DataFrame
# because it does not have a geometry column anymore. Also note that
# it is assumed the `'geom'` column is correctly datatyped.
gdf['geom'] = gdf['geometry'].apply(wkb_hexer)

# Delete extra WKT geometry colum.
del gdf['geometry']

# Empty column.
del gdf['WIDTH']

table_name = 'sj_sidewalks_poly''

# Connect to database using a context manager
gdf.to_sql(table_name, con=conn, if_exists='append', index=False)


#C onvert the geom column to a geometry.
sql = """ALTER TABLE %s
ALTER COLUMN geom TYPE Geometry(POLYGON, 4326)
                  USING ST_SetSRID(geom::Geometry, 4326)""" % (table_name)
conn.execute(sql)
# -

# ### Make sure PostGIS is running right
# You may need to do a series of these commands to get your postgis running the right SFCGAL pacakges.
#
# ```CREATE extension postgis_sfcgal```
#
# ```SET postgis.backend = sfcgal;```
#
# ```SET postgis.backend = geos;```
#
# Run this command to test to see if SFCGAL is working:
# ```SELECT ST_ApproximateMedialAxis(ST_GeomFromText('POLYGON (( 190 190, 10 190, 10 10, 190 10, 190 20, 160 30, 60 30, 60 130, 190 140, 190 190 ))'));```

# ## PostGIS SQL Processing Commands  

# Take the 44k sidewalk polygons and merge them all into one big geometry
# There are lots of sidewalk polygons next to each other that refer to connected sidewalks.  
# During this step we lose the ADACOMPLY, COVERED, and other fields.
sql = """CREATE TABLE union_sideys as SELECT ST_UNION(geom) as geom FROM sj_sidewalks_poly;"""
conn.execute(sql)

# Break out the single polygon into the individual polygons.  This produces about 7375 total rows.
sql = """CREATE TABLE unjoined_sideys as SELECT (ST_Dump(geom)).geom AS geom FROM union_sideys;"""
conn.execute(sql)

#Give the table some ids because of the next step.
#https://dba.stackexchange.com/questions/20801/most-efficient-way-to-add-a-serial-column-to-a-huge-table
sql = """ALTER TABLE unjoined_sideys ADD column id bigserial;"""
conn.execute(sql)

# Create the empty table that will contain the sidewalk linestrings.
sql = """CREATE TABLE approx_medial (id bigserial, LINESTRING geometry); """
conn.execute(sql)

error_rows = []

# +
# Individually process each of the 7375 total polygons.  This algorithm mostly produces the
# middle line of the polygons represented by a set of lines which is close to what we want.
# https://postgis.net/docs/manual-2.2/ST_ApproximateMedialAxis.html

# Unfortunately 67 of the 7.3k rows will fail this algorithm since: 
# "ERROR: straight skeleton of Polygon with touching interior rings is not implemented"
# This method takes a while, maybe around 30 minutes.  It can probably be optimized further.

for row in range(0,7376):
    sql = """INSERT INTO approx_medial (SELECT id, ST_ApproximateMedialAxis(geom) AS linestring FROM unjoined_sideys WHERE id = %s);""" % (row)
    try:
        conn.execute(sql)
    except:
        error_rows.append(row)
# -

error_rows = [1165, 1166, 1422, 1680, 1696, 2024, 2025, 2247, 2381, 2382, 2874, 2979, 3255, 3256, 3537, 3538, 3539, 3540, 3541, 3542, 3543, 3544, 3545, 3546, 3547, 3548, 3549, 3868, 3869, 3941, 3942, 4120, 4121, 4465, 4630, 4631, 4632, 4633, 4634, 4753, 4754, 4956, 4957, 4958, 4959, 4960, 5182, 5183, 5184, 5185, 5790, 5791, 6043, 6173, 6174, 6175, 6229, 6471, 6472, 6542, 6613, 6644, 6703, 6774, 6866, 6867, 7046]
# We'll save the broken geometrys if we want a reference for us to do manually.


# The ApproximateMedialAxis lines created point in different directions,
# this step simplifies them.
sql = '''CREATE TABLE merged_approx_medial as SELECT id, ST_LineMerge(linestring) AS linestring FROM approx_medial;'''
conn.execute(sql)

# This specific simplify value (.00001) works perfectly for this case. Determined through
# much testing.
sql = """CREATE TABLE approx_medial_simplified AS 
SELECT ST_Simplify(linestring,.00001) as linestring, id 
FROM merged_approx_medial"""
conn.execute(sql)

#
# ### Get rid of barbs
# At this point in the process we have pretty good backbones of the sidewalk network but there are these little barbs poking out.  We should remove as many as possible to save us work during the import process

# +
# Get the length of all the lilnestrings and then only keep the linestrings
# that are longer than .000038.  This number was found through careful testing

sql = """CREATE TABLE sim_length as 
SELECT (ST_Dump(linestring)).geom AS linestring, 
ST_Length((ST_Dump(linestring)).geom) AS length, ID 
FROM approx_medial_simplified;"""
conn.execute(sql)

sql = """CREATE TABLE merged_again AS 
SELECT ST_LineMerge(l.linestring) AS linestring 
FROM (SELECT ST_Union(k.linestring) AS linestring 
FROM (SELECT * FROM sim_length WHERE length > 0.000038) k) l;"""
conn.execute(sql)
# -

sql = """CREATE TABLE mtwice as SELECT (ST_Dump(ST_LineMerge(ST_Union(linestring)))).geom as linestring FROM merged_processed_again;"""
conn.execute(sql)

# ### Deal with tiny holes
# At this point the barbs are gone but also these tiny holes appeared in the sidewalks where the barbs were.  For almost all the data we can just connect the end of one line to the beginning of another if a new line is within .000038.
#

#Code creates starts and ends table, the start of lines and end of lines saved as points.
sql = """CREATE TABLE starts AS SELECT ST_StartPoint(linestring) AS point FROM mtwice;"""
conn.execute(sql)
sql = """CREATE TABLE ends AS SELECT ST_EndPoint(linestring) AS point FROM mtwice;"""
conn.execute(sql)

# +
# Find the distances between the start points and the nearest end points,
# start points and the nearest other start points,
# end points and the nearest other end points.
sql = """CREATE TABLE start_to_end_short AS SELECT DISTINCT ON (starts.point) starts.point AS sp, ends.point AS ep, ST_Distance(starts.point, ends.point) AS dist
FROM starts AS starts, ends AS ends WHERE ST_DWithin(starts.point, ends.point, 0.000038);"""
conn.execute(sql)

sql = """CREATE TABLE start_to_start_short AS SELECT DISTINCT ON (a.point) a.point AS sp, b.point AS ep, ST_Distance(a.point, b.point) AS dist
FROM starts a, starts b WHERE ST_DWithin(a.point, b.point, 0.000038);"""
conn.execute(sql)

sql = """CREATE TABLE end_to_end_short AS SELECT DISTINCT ON (a.point) a.point AS sp, b.point AS ep, ST_Distance(a.point, b.point) AS dist
FROM ends a, ends b WHERE ST_DWithin(a.point, b.point, 0.000038);"""
conn.execute(sql)

# +
# Insert the mini-lines to fill in the gaps.  Only take distances greater than 0.
sql = """INSERT INTO mtwice (linestring) SELECT ST_MakeLine(sp, ep) AS linestring FROM start_to_end_short WHERE dist != 0;"""
conn.execute(sql)

sql = """INSERT INTO mtwice (linestring) SELECT ST_MakeLine(sp, ep) AS linestring FROM start_to_start_short WHERE dist != 0;"""
conn.execute(sql)

sql = """INSERT INTO mtwice (linestring) SELECT ST_MakeLine(sp, ep) AS linestring FROM end_to_end_short WHERE dist != 0;"""
conn.execute(sql)
# -

# Merge and simplify the mini lines to the rest of the lines one final time.
sql = """CREATE TABLE mthrice as 
SELECT ST_Simplify((ST_Dump(ST_LineMerge(ST_Union(linestring)))).geom,.00001) as linestring FROM mtwice;"""
conn.execute(sql)

# +
sql = """SELECT * FROM mthrice"""
sw_m = gpd.read_postgis(sql, conn, geom_col='linestring')
sw_m.crs = {'init': 'epsg:4326'}


# Properly tag the sidewalk data
sw_m['surface'] = 'paved'
sw_m['highway'] = 'footway'
sw_m['footway'] = 'sidewalk'

# Get the centroids of the lines
sw_m['geometry'] = sw_m.centroid
sw_m = sw_m.set_geometry('geometry')

# +
# Pull little census blocks of San Jose using VTA TAZs.
tazs = gpd.read_file('vta_method/vta_method.shp')
tazs = tazs.to_crs({'init': 'epsg:4326'})

tazs = tazs.loc[tazs['geometry'].geom_type == 'Polygon',]
tazs = tazs[['VTA_TAZ','geometry']]
tazs = tazs.set_geometry('geometry')
# -

# Group the linestrings by doing a spatial join with the centroids
# thus each linestring has a zone tag for the right area.
joined = gpd.sjoin(tazs,sw_m, how="right", op='contains')
del joined['index_left']
del joined['geometry']
joined = joined.set_geometry('linestring')
zones = joined['VTA_TAZ'].unique()

# Query and create separate geojsons for each zone of sidewalks
for zone in zones:
    filename = 'sidewalk_zones/sw_zone_%d' % zone + '.geojson'
    z_query = "VTA_TAZ==%d" % (zone)
    joined.query(z_query)[['linestring','surface','highway','footway']].to_file(filename, driver='GeoJSON')

# ## Get and use osmizer
#
# When you import the geojson into to get rid of duplicate nodes
#
# https://github.com/OpenSidewalks/osmizer
#
# Run this command using osmizer to produce final .osm files:
# ```find . -name \*.geojson -exec sh -c 'osmizer convert sidewalks {} $(basename {} .geojson).osm' \;```
#
# ```rm -rf *.geojson```

# +
# Only output the tazs that have sidewalks within them.  
# These will become our tasks.
# Import URL helps with JOSM

tazs = tazs.loc[tazs['VTA_TAZ'].isin(zones),]

def formatter(word):
    initial_url = 'https://codeforsanjose.github.io/OSM-SouthBay/SJ_Sidewalks/import_data/'
    return initial_url + 'sw_zone_' + word + '.osm'

tazs['import_url']= tazs['VTA_TAZ'].astype(str).apply(formatter)

tazs.to_file('sidewalks_tazs.geojson', driver='GeoJSON')
tazs.head()

StarcoderdataJupyterScriptsDedupFiltered