Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Highway length to building count ratio for NYC
WITH country AS (
SELECT
osm_layers.name as osm_name,
osm_layers.all_tags AS osm_tags,
(SELECT tags.value FROM UNNEST(all_tags) as tags WHERE tags.key = 'admin_level') as admin_level,
osm_layers.geometry AS geometry
FROM `openstreetmap-public-data-dev.osm_planet.osm_layers` AS osm_layers
WHERE osm_layers.name='boundary-administrative'
AND EXISTS(SELECT tags.value FROM UNNEST(all_tags) as tags WHERE tags.key = 'name' and tags.value = 'United States')
),
country_city AS (
SELECT
osm_layers.name as osm_name,
osm_layers.all_tags AS osm_tags,
(SELECT tags.value FROM UNNEST(all_tags) as tags WHERE tags.key = 'admin_level') as admin_level,
osm_layers.geometry AS geometry
FROM `openstreetmap-public-data-dev.osm_planet.osm_layers` AS osm_layers, country
WHERE osm_layers.name='boundary-administrative'
AND EXISTS(SELECT tags.value FROM UNNEST(all_tags) as tags WHERE tags.key = 'name' and tags.value='New York')
AND EXISTS(SELECT tags.value FROM UNNEST(all_tags) as tags WHERE tags.key = 'place' and tags.value='city')
AND ST_DWITHIN(country.geometry, osm_layers.geometry, 0)
),
country_city_buildings AS (
SELECT
osm_layers.*
FROM `openstreetmap-public-data-dev.osm_planet.osm_layers` AS osm_layers, country_city
WHERE osm_layers.name='building'
AND ST_DWITHIN(country_city.geometry, osm_layers.geometry, 0)
),
country_city_highways AS (
SELECT
osm_layers.*
FROM `openstreetmap-public-data-dev.osm_planet.osm_layers` AS osm_layers, country_city
WHERE osm_layers.name='highway'
AND ST_DWITHIN(country_city.geometry, osm_layers.geometry, 0)
),
country_city_buildings_agg AS (
SELECT
count(1) as count,
ST_GEOHASH(ST_SNAPTOGRID(ST_CENTROID(geometry), 0.01)) AS geohash
FROM country_city_buildings
GROUP BY 2
),
country_city_buildings_agg_points AS (
SELECT
count,
ST_GEOGPOINTFROMGEOHASH(geohash) AS geometry
FROM country_city_buildings_agg
),
country_city_buildings_agg_cells AS (
SELECT
count,
ST_GEOHASH(geometry) as geohash,
(SELECT
ST_MAKEPOLYGON(ST_MAKELINE(ARRAY_AGG(geom)))
FROM UNNEST(ARRAY[
ST_GEOGPOINT(ST_X(geometry)-0.25/50, ST_Y(geometry)-0.25/50),
ST_GEOGPOINT(ST_X(geometry)-0.25/50, ST_Y(geometry)+0.25/50),
ST_GEOGPOINT(ST_X(geometry)+0.25/50, ST_Y(geometry)+0.25/50),
ST_GEOGPOINT(ST_X(geometry)+0.25/50, ST_Y(geometry)-0.25/50)
]) as geom
) as geometry
FROM country_city_buildings_agg_points
),
buildings_area AS (
SELECT
geohash,
count(1) AS count,
SUM(ST_AREA(ST_INTERSECTION(buildings.geometry,cells.geometry))) AS area
FROM country_city_buildings_agg_cells AS cells, country_city_buildings AS buildings
WHERE ST_INTERSECTS(buildings.geometry,cells.geometry)
GROUP BY 1
),
highways_length AS (
SELECT
geohash,
count(1) AS count,
SUM(ST_LENGTH(ST_INTERSECTION(highways.geometry,cells.geometry))) AS length
FROM country_city_buildings_agg_cells AS cells, country_city_highways AS highways
WHERE ST_INTERSECTS(highways.geometry,cells.geometry)
GROUP BY 1
)
SELECT
cells.geometry,
buildings_area.count AS buildings_count,
buildings_area.area AS buildings_area,
highways_length.length AS highways_length,
highways_length.count AS highways_count,
(CASE WHEN highways_length.length=0 THEN NULL ELSE highways_length.length END)
/
(CASE WHEN buildings_area.area=0 THEN NULL ELSE buildings_area.area END)
AS highways2buildings
FROM buildings_area, highways_length, country_city_buildings_agg_cells as cells
WHERE buildings_area.geohash = highways_length.geohash
AND buildings_area.geohash = cells.geohash;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.