Highway length to building count ratio for NYC

OSM-NYC-highways-to-buildings.sql

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;