Geo-Rails Part 8: ZCTA Lookup, A Worked Example

actual

% rails generate model Zcta zcta:integer region:polygon
class CreateZctas < ActiveRecord::Migration
  def change
    create_table :zctas do |t|
      t.integer :zcta
      t.polygon :region
      t.timestamps
    end
  end
end
class CreateZctas < ActiveRecord::Migration
  def change
    create_table :zctas do |t|
      t.integer :zcta
      t.polygon :region, :srid => 3785
    end
    change_table :zctas do |t|
      t.index :region, :spatial => true
    end
  end
end
% rake db:migrate
class Zcta < ActiveRecord::Base
 FACTORY = RGeo::Geographic.simple_mercator_factory
 set_rgeo_factory_for_column(:region, FACTORY.projection_factory)
end
require 'rgeo-shapefile'
RGeo::Shapefile::Reader.open('tl_2010_53_zcta510.shp',
    :factory => Zcta::FACTORY) do |file|
  file.each do |record|
    zcta = record['ZCTA5CE10'].to_i
    # The record geometry is a MultiPolygon. Iterate
    # over its parts.
    record.geometry.projection.each do |poly|
      Zcta.create(:zcta => zcta, :region => poly)
    end
  end
end
class Zcta < ActiveRecord::Base

  # ...

  EWKB = RGeo::WKRep::WKBGenerator.new(:type_format => :ewkb,
    :emit_ewkb_srid => true, :hex_format => true)

  def self.containing_latlon(lat, lon)
    ewkb = EWKB.generate(FACTORY.point(lon, lat).projection)
    where("ST_Intersects(region, ST_GeomFromEWKB(E'\\\\x#{ewkb}'))")
  end

end
class Zcta < ActiveRecord::Base

  # ...

  def self.containing_geom(geom)
    ewkb = EWKB.generate(FACTORY.project(geom))
    where("ST_Intersects(region, ST_GeomFromEWKB(E'\\\\x#{ewkb)}'))")
  end

end
class ZctaController

  def lookup
    lat = params[:lat].to_f
    lon = params[:lon].to_f
    zcta = Zcta.containing_latlon(lat, lon).first
    render(:json => {:lat => lat, :lon => lon,
      :zcta => zcta ? zcta.zcta : nil})
  end

end
MAX_SIZE = 500
MAX_DEPTH = 12

require 'rgeo-shapefile'

# Handle a geometry of any type
def handle_geometry(depth, geom, zcta)
  case geom
  when ::RGeo::Feature::Polygon
    handle_polygon(depth, geom, zcta)
  when ::RGeo::Feature::MultiPolygon
    geom.each{ |polygon| handle_polygon(depth, polygon, zcta) }
  end
end

# Handle a polygon
def handle_polygon(depth, polygon, zcta)
  # Check the number of sides. We'll combine the number of sides for
  # the "outer edge" and any "holes" that the polygon might have.
  # A polygon boundary consists of a LineString that is closed, so
  # the first and last points are the same. Therefore, to count the
  # sides, count the number of vertices and subtract 1.
  sides = polygon.exterior_ring.num_points - 1
  polygon.interior_rings.each{ |ring| sides += ring.num_points - 1 }
  if depth >= MAX_DEPTH || sides <= MAX_SIZE
    # The polygon is small enough, or we recursed as far as we're
    # willing. Just add the polygon.
    Zcta.create(:zcta => zcta, :region => polygon)
  else
    # Split the polygon 4-to-1 and recurse
    depth = depth + 1
    # Find the bounding box for the polygon
    envelope = polygon.envelope.exterior_ring
    p1 = envelope.point_n(0)
    p2 = envelope.point_n(2)
    min_x = p1.x
    max_x = p2.x
    min_x, max_x = max_x, min_x if min_x > max_x
    min_y = p1.y
    max_y = p2.y
    min_y, max_y = max_y, min_y if min_y > max_y
    # dx and dy are the size of the bounding box.
    # cx and cy are the center point.
    dx = max_x - min_x
    dy = max_y - min_y
    cx = (min_x + max_x) * 0.5
    cy = (min_y + max_y) * 0.5
    # Check the aspect ratio of the bounding box. If it's very wide
    # or very tall, then only split in half. Otherwise, split in 4.
    if dy > dx * 2
      # The bounding box is tall, so split in half
      handle_quadrant(depth, polygon, min_x, min_y, max_x, cy, zcta)
      handle_quadrant(depth, polygon, min_x, cy, max_x, max_y, zcta)
    elsif dx > dy * 2
      # The bounding box is wide, so split in half
      handle_quadrant(depth, polygon, min_x, min_y, cx, max_y, zcta)
      handle_quadrant(depth, polygon, cx, min_y, max_x, max_y, zcta)
    else
      # The bounding box is close to square so split in four
      handle_quadrant(depth, polygon, min_x, min_y, cx, cy, zcta)
      handle_quadrant(depth, polygon, cx, min_y, max_x, cy, zcta)
      handle_quadrant(depth, polygon, min_x, cy, cx, max_y, zcta)
      handle_quadrant(depth, polygon, cx, cy, max_x, max_y, zcta)
    end
  end
end

# Take a polygon and a box. Run the algorithm on the part of the
# polygon that falls within the box.
def handle_quadrant(depth, polygon, min_x, min_y, max_x, max_y, zcta)
  # We do this by creating a rectangle for the box, and computing
  # the intersection with the input polygon. The result could be a
  # polygon, a MultiPolygon, or an empty geometry.
  factory = Zcta::FACTORY.projection_factory
  box = factory.polygon(factory.linear_ring([
    factory.point(min_x, min_y),
    factory.point(min_x, max_y),
    factory.point(max_x, max_y),
    factory.point(max_x, min_y)]))
  handle_geometry(depth, polygon.intersection(box), zcta)
end

# The main shapefile reader.
RGeo::Shapefile::Reader.open('tl_2010_53_zcta510.shp',
    :factory => Zcta::FACTORY) do |file|
  file.each do |record|
    # For each MultiPolygon, analyze it and add to the database
    handle_geometry(0, record.geometry.projection,
      record['ZCTA5CE10'].to_i)
  end
end

Geo-Rails Part 8: ZCTA Lookup, A Worked Example

ZCT-what?

Setting up our database to hold ZCTA data

Obtaining ZCTA data and populating our database

Running ZCTA queries

Segmenting polygons for improved performance

Where to go from here

Dialogue & Discussion