Skip to content

Instantly share code, notes, and snippets.

@andrewxhill
Created October 9, 2012 15:17
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save andrewxhill/3859464 to your computer and use it in GitHub Desktop.
Save andrewxhill/3859464 to your computer and use it in GitHub Desktop.
Create an equal area grid over a buffer of points
-- 1) Create the desired buffer of 50km around your points, buffs
-- 2) Create a reprojected envelope around your buffs, area
-- 3) Create a rectangular grid inside the envelope, grid
-- 4) Select the centroid of each gridcell and reproject it as a a WGS84 point
-- The equal area Lambert projection for North America is documented here,
-- http://spatialreference.org/ref/sr-org/7314/
-- and is stored in your spatial ref system as, 97314
WITH buffs AS ( --creates a resource called 'area' that will act like any existing table
SELECT
ST_Buffer(
the_geom::geography,
50000)::geometry AS buffers,
name
FROM uofa_4km_nad83
),
area AS ( -- same as buffs, create a second resource
SELECT
ST_Envelope( -- returns a polygon geometry that covers everything included
ST_Transform( -- changes the porojection
ST_Union(buffers), -- joins all returned geometries into one geom
97314 -- the equal-area Lamber projection id for ST_Transform
)
) as envelope -- call this column, 'the_geom'
FROM buffs --using buffs from above as our table
),
grid AS (
SELECT
CDB_RectangleGrid( -- A cartodb function to return rect grid polygons
envelope, --It will return the grid in whatever projection the input geometry is in
1000, --width
1000 --height
) cell --call the column 'cell'
FROM area
)
SELECT
ST_Centroid( --Returns the centroid of any geometry as a point
ST_Transform( -- transform our equal area projection back into WGS84
cell,
4326
)
) as the_geom
FROM grid
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment