Using PostGIS

creating template_postgis

http://geospatial.nomad-labs.com/2006/12/24/postgis-template-database/

sliver problems - part 1

This idiom can help filter out slivers based on size.
This query selects the features that intersect,
it also creates an intersection of those features,
and returns the resultset.

This can be used to detect and filter out slivers.

select
	zd.zoning_sim as zone,
	st_area(st_intersection(cl.the_geom, zd.the_geom)) as area
from
	v_citylots cl,
	planning.zoning_districts zd
where cl.blklot = '3514003'
and st_intersects(cl.the_geom, zd.the_geom)
order by st_area(st_intersection(cl.the_geom, zd.the_geom)) desc

sliver problems - part 2

If a simple area threshhold fails to allow sliver filtering
we might try a metric from ecology, the area to perimeter ratio.
Here is a query:

select
	zd.zoning_sim as zone,
	st_area(st_intersection(cl.the_geom, zd.the_geom)) as area,
	st_length(st_boundary(st_intersection(cl.the_geom, zd.the_geom))) as length,
	st_area(st_intersection(cl.the_geom, zd.the_geom)) / st_length(st_boundary(st_intersection(cl.the_geom, zd.the_geom))) as ap_ratio,
	st_length(st_boundary(st_intersection(cl.the_geom, zd.the_geom))) / st_area(st_intersection(cl.the_geom, zd.the_geom)) as pa_ratio
from
	v_citylots cl,
	planning.zoning_districts zd
where cl.blklot = '3514003'
and st_intersects(cl.the_geom, zd.the_geom)
order by st_area(st_intersection(cl.the_geom, zd.the_geom)) desc

And here are the results:

zone	area	length	a2p	p2a
C-3-G	1756.165	172.792	10.163	0.098
NCT-3	0.618	66.608	0.009	107.696

The first row is a fairly normal city lot.
The second row is the sliver.
As you can see, the difference is an order of magnitude which is what we would want for a threshhold.

convert multi geometry to single geometry

These statements create a single geometry table from a multi geometry table.
The analogous arcmap operation is called explode (I think).
This is useful to help improve the rendering performance in WMS.

-- Create the single geometry is from the multi geometry.
create table my_single_geom as
select gid, objectid, distance, area, length, label, (ST_Dump(the_geom)).geom as the_geom
from my_multi_geom

-- Insert a row into geometry_columns.
-- Without this, its not a postgis layer.
-- A single geometry polygon is of type is POLYGON.
-- A multiple geometry polygon is of type is MULTIPOLYGON.
INSERT INTO geometry_columns (
    f_table_catalog,
    f_table_schema,
    f_table_name character,
    f_geometry_column,
    coord_dimension,
    srid,
    type)
VALUES (
    '',
    '',
    'my_single_geom'
    'the_geom',
    2,
    900913,
    'POLYGON'
)