PostGIS: Up and Running

Regina Obe and Leo Hsu



Buy our books! at http://www.postgis.us



PGOpen promo codes for 44% off PostGIS In Action 2nd Edition postgrescft



Tutorial material at http://www.postgis.us/pgopen2014

Agenda

Preliminaries

Geometry

Spatial Reference Systems

Geography

Geocoding

3D Geometries


Tantalizing SQL Constructs

Importing Data

Tools of the Trade

Geoid

Ellipsoid

OS2PGSQL and other OpenStreetMap loaders

  • Osm2pgsql http://wiki.openstreetmap.org/wiki/Osm2pgsql is one of several tools not packaged with PostGIS but useful for loading OSM data into PostGIS format. Supported on all platforms and was one of the first

  • Imposm http://imposm.org - another popular for Linux/Mac OSX - no binaries for windows. Many say it is faster than osm2pgsql.

  • ogr2ogr Part of GDAL toolkit. Binaries available for many platforms and often already installed if you installed PostGIS on Linux.

3D Surfaces

We are using PostGIS ST_AsX3D function combined with http://www.x3dom.org to render some of these in HTML

Warning 3D stuff is still not very tested yet.

Closed Polyhedral Surface

SELECT ST_GeomFromText('POLYHEDRALSURFACE Z( 
   ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0))
 , ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0))
 , ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)) 
 , ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)) 
 , ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0))
 , ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )');

Closed Polyhedral Surface In X3D

SELECT ST_AsX3D('POLYHEDRALSURFACE Z( 
   ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0))
 , ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0))
 , ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)) 
 , ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)) 
 , ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0))
 , ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )'::geometry);

Output Closed PolyhedralSurface with X3D

SELECT '
 ' || 
ST_AsX3D('POLYHEDRALSURFACE Z( 
   ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0))
 , ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0))
 , ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)) 
 , ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)) 
 , ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0))
 , ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )'::geometry) || '';

Closed TIN

SELECT 'TIN(
            ((0 0 0, 1 0 0, 0 1 0, 0 0 0)),
            ((0 0 0, 1 0 0, 0 0 1, 0 0 0)),
            ((0 0 0, 0 0 1, 0 1 0, 0 0 0)),
            ((0 0 1, 1 0 0, 0 1 0, 0 0 1))
            )'::geometry;

Closed TIN with X3D

SELECT '
  ' 
    || replace(ST_AsX3D('TIN(
    ((0 0 0, 1 0 0, 0 1 0, 0 0 0)),
    ((0 0 0, 1 0 0, 0 0 1, 0 0 0)),
    ((0 0 0, 0 0 1, 0 1 0, 0 0 0)),
    ((0 0 1, 1 0 0, 0 1 0, 0 0 1))
            )'::geometry),  
   '<IndexedTriangleSet', '<IndexedTriangleSet solid="false" ') || '';

Road network in X3D

SELECT string_agg('
 
' || 
ST_AsX3D(ST_Force_3D(geom)) || '', '') As x3d
FROM (SELECT row_number() over(ORDER BY ST_Length(way) DESC) As rn,
 name, count(*) OVER() As cnt, way As geom
 FROM po.planet_osm_line 
   WHERE ST_Intersects(ST_Transform(
    ST_Buffer(
       ST_Point(-87.627,41.8819)::geography,
       100)::geometry,900913
        ) , way)  ) As f;

SFCGAL ST_Extrude convert 2D Buildings to 3D

SELECT string_agg('
 
' || 
ST_AsX3D(ST_Extrude(geom,0,0,50)) || '', '')  As x3d
FROM (SELECT row_number() over(ORDER BY ST_Area(p.way) DESC) As rn,
 name, count(*) OVER() As cnt, p.way As geom
 FROM po.planet_osm_polygon AS p 
    INNER JOIN (SELECT way FROM po.planet_osm_polygon
 WHERE name ILIKE 'Trump%') As ref ON ST_DWithin(ref.way, p.way,150)
 WHERE p.building > ''
  ) As f;

SFCGAL ST_Extrude Output

Tiger Geocoder: Loading data

Create copy of a template record

Here I use windows. If on Linux or mac use sh

INSERT INTO tiger.loader_platform(os, declare_sect, pgbin, 
	wget, unzip_command, psql, path_sep, 
		   loader, environ_set_command, county_process_command)
SELECT 'mydb', declare_sect, pgbin, wget, unzip_command, psql, path_sep, 
	   loader, environ_set_command, county_process_command
  FROM tiger.loader_platform
  WHERE os = 'windows';

Use your favorite editor to edit the declare_sect.

Generate Nation files load script

SELECT Loader_Generate_Nation_Script('mydb');

Then run the generated script from command line.

Tip 1: to reduce load skip blockgroups

If you don't need tabblocks, block groups, and tracts for statistics, disable load. These aren't used by core geocoder functionality. Disable the load.

UPDATE tiger.loader_lookuptables 
	SET load = false WHERE table_name IN('bg', 'tabblock', 'tract');

Generate State load script

SELECT Loader_Generate_Script('{IL}'::text[], 'mydb');

Tip 2: Don't need whole state edit down to county

If you don't need a whole state of data, edit the generated script.

For example faces, edges, featnames -- change 17 to 17031 (if just want Cook county)

Then run the generated script from command line.

Install missing indexes and analyze

SELECT Install_Missing_Indexes();

don't forget to update stats

analyze verbose;

Geocoding

SELECT pprint_addy(addy), 
 ST_AsText(ST_SnapToGrid(geomout,0.0001)), rating 
FROM geocode('333 North Dearborn Street, Chicago, IL 60654',1);
             pprint_addy              |        st_astext        | rating
--------------------------------------+-------------------------+--------
 333 N Dearborn St, Chicago, IL 60654 | POINT(-87.6294 41.8885) |      0

Geocoding Intersections

Intersection Name on Lake St

SELECT pprint_addy(addy), 
 ST_AsText(ST_SnapToGrid(geomout,0.0001)) As loc, rating 
FROM geocode_intersection('Lake', 'Clark', 'IL', 'Chicago');
           pprint_addy            |          loc           | rating
----------------------------------+------------------------+--------
 101 W Lake St, Chicago, IL 60601 | POINT(-87.631 41.8857) |      5
 100 W Lake St, Chicago, IL 60661 | POINT(-87.631 41.8857) |      5
 W Lake St, Chicago, IL           | POINT(-87.631 41.8857) |      5

Geocoding Intersections

Intersection Name on Clark St

SELECT pprint_addy(addy), 
 ST_AsText(ST_SnapToGrid(geomout,0.0001)) As loc, rating 
FROM geocode_intersection('Clark', 'Lake', 'IL', 'Chicago');
            pprint_addy            |          loc           | rating
-----------------------------------+------------------------+--------
 200 N Clark St, Chicago, IL 60601 | POINT(-87.631 41.8857) |      5
 201 N Clark St, Chicago, IL 60601 | POINT(-87.631 41.8857) |      5
 198 N Clark St, Chicago, IL 60601 | POINT(-87.631 41.8857) |      5
 199 N Clark St, Chicago, IL 60601 | POINT(-87.631 41.8857) |      5

Geocoding Intersections, include zip and num matches

You get better performance if you provide Zip and max returns

SELECT pprint_addy(addy), 
 ST_AsText(ST_SnapToGrid(geomout,0.00001)) As loc, rating 
FROM geocode_intersection('Clark', 'Lake', 'IL', 'Chicago', '60601', 1);
            pprint_addy            |            loc            | rating
-----------------------------------+---------------------------+--------
 200 N Clark St, Chicago, IL 60601 | POINT(-87.63098 41.88574) |      5

Reverse Geocoding

addy is an array of possible normalized addresses of a point location. You can expand aliases with unnest.

SELECT pprint_addy(unnest(rc.addy))
 FROM reverse_geocode(ST_Point(-87.63098, 41.88574) ) As rc;

If you are running PostGIS < 2.1.4, you will be missing the directional N.

           pprint_addy
---------------------------------
 102 N Lake St, Chicago, IL 60661
 200 N Clark St, Chicago, IL 60601

Links of Interest

Promo Codes

44% off PostGIS In Action 2nd Edition postgrescft