Simple SQL GIS
31 Jul 2019And, late on a Friday afternoon, the plaintive cry was heard!
Anyone got a KML/Shapefile of B.C. elxn boundaries that follows the water (Elections BC's KML has ridings going out into the sea)
— Chad Skelton (@chadskelton) November 16, 2012
And indeed, into the sea they do go!
And ‘lo, the SQL faeries were curious, and gave it a shot!
##### Commandline OSX/Linux #####
# Get the Shape files
# http://www.elections.bc.ca/index.php/voting/electoral-maps-profiles/
wget http://www.elections.bc.ca/docs/map/redis08/GIS/ED_Province.exe
# Exe? No prob, it's actually a self-extracting ZIP
unzip ED_Province
# Get a PostGIS database ready for the data
createdb ed_clip
psql -c "create extension postgis" -d ed_clip
# Load into PostGIS
# The .prj says it is "Canada Albers Equal Area", but they
# lie! It's actually BC Albers, EPSG:3005
shp2pgsql -s 3005 -i -I ED_Province ed | psql -d ed_clip
# We need some ocean! Use Natural Earth...
# http://www.naturalearthdata.com/downloads/
wget http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_ocean.zip
unzip ne_10m_ocean.zip
# Load the Ocean into PostGIS!
shp2pgsql -s 4326 -i -I ne_10m_ocean ocean | psql -d ed_clip
# OK, now we connect to PostGIS and start working in SQL
psql -e ed_clip
-- How big is the Ocean table?
SELECT Count(*) FROM ocean;
-- Oh, only 1 polygon. Well, that makes it easy...
-- For each electoral district, we want to difference away the ocean.
-- The ocean is a one big polygon, this will take a while (if we
-- were being more subtle, we'd first clip the ocean down to
-- a reasonable area around BC.)
CREATE TABLE ed_clipped AS
SELECT
CASE
WHEN ST_Intersects(o.geom, ST_Transform(e.geom,4326))
THEN ST_Difference(ST_Transform(e.geom,4326), o.geom)
ELSE ST_Transform(e.geom,4326)
END AS geom,
e.edabbr,
e.edname
FROM ed e, ocean o;
-- Check our geometry types...
SELECT DISTINCT ST_GeometryType(geom) FROM ed_clipped;
-- Oh, they are heterogeneous. Let's force them all multi
UPDATE ed_clipped SET geom = ST_Multi(geom);
# Dump the result out of the database back into shapes
pgsql2shp -f ed2009_ocean ed_clip ed_clipped
zip ed2009_ocean.zip ed2009_ocean.*
mv ed2009_ocean.zip ~/Dropbox/Public/
No more districts in oceans!
And the faeries were happy, and uploaded their polygons!
Update: And the lamentations ended, and the faeries also rejoiced.
@pwramsey OK, that's frickin' amazing! Thank you! Thought I was in store for hours spent editing polygons by hand in Google Earth.
— Chad Skelton (@chadskelton) November 17, 2012