Importing Data

CSV Import: Create table

This needs to mirror structure of your import file

CREATE TABLE data.london_stations
(
   station character varying(150), 
   os_x integer, 
   os_y integer, 
   latitude numeric(10,6),
   longitude numeric(10,6),
   zone varchar(20), postcode varchar(50),
   constraint pk_long_stations primary key(station) 
) ;

CSV Import: Load Data

For this connect to psql -- either via pgAdmin plugin or from OS console

This step already done for you if connecting via pgAdmin psql

\connect workshop_1

Load data

\cd path_to_folder
\copy data.london_stations FROM 'London stations.csv' DELIMITER ',' CSV HEADER;

shp2pgsql: Set environment variables from commandline

From commandline set useful variables -- use SET on windows export on Linux. Setting path is usually not necessary for Linux/MacOSX

SET PATH=%PATH%;c:/Program Files/PostgreSQL/9.4/bin
SET PGUSER=postgres
SET PGPASSWORD=whatever
SET PGHOST=localhost
SET PGPORT=5432
SET PGDATABASE=workshop_1

shp2pgsql: get list of commands

shp2pgsql

RELEASE: 2.2.0dev (r13180)
USAGE: shp2pgsql [<options>] <shapefile> [[<schema>.]<table>]
OPTIONS:
  -s [<from>:]<srid> Set the SRID field. Defaults to 0.
      Optionally reprojects from given SRID (cannot be used with -D).
 (-d|a|c|p) These are mutually exclusive options:
     -d  Drops the table, then recreates it and populates
         it with current shape file data.
     -a  Appends shape file into current table, must be
         exactly the same table schema.
     -c  Creates a new table and populates it, this is the
         default if you do not specify any options.
     -p  Prepare mode, only creates the table.
  -g <geocolumn> Specify the name of the geometry/geography column
      (mostly useful in append mode).
  -D  Use postgresql dump format (defaults to SQL insert statements).
  -e  Execute each statement individually, do not use a transaction.
      Not compatible with -D.
  -G  Use geography type (requires lon/lat data or -s to reproject).
  -k  Keep postgresql identifiers case.
  -i  Use int4 type for all integer dbf fields.
  -I  Create a spatial index on the geocolumn.
  -m <filename>  Specify a file containing a set of mappings of (long) column
     names to 10 character DBF column names. The content of the file is one or
     more lines of two names separated by white space and no trailing or
     leading space. For example:
         COLUMNNAME DBFFIELD1
         AVERYLONGCOLUMNNAME DBFFIELD2
  -S  Generate simple geometries instead of MULTI geometries.
  -t <dimensionality> Force geometry to be one of '2D', '3DZ', '3DM', or '4D'
  -w  Output WKT instead of WKB.  Note that this can result in
      coordinate drift.
  -W <encoding> Specify the character encoding of Shape's
      attribute column. (default: "UTF-8")
  -N <policy> NULL geometries handling policy (insert*,skip,abort).
  -n  Only import DBF file.
  -T <tablespace> Specify the tablespace for the new table.
      Note that indexes will still use the default tablespace unless the
      -X flag is also used.
  -X <tablespace> Specify the tablespace for the table's indexes.
      This applies to the primary key, and the spatial index if
      the -I flag is used.
  -?  Display this help screen.

Guessing at spatial projection from building_footprint.prj file

PROJCS["NAD_1983_StatePlane_California_III_FIPS_0403_Feet",
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.2572221]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199]],PROJECTION["Lambert_Conformal_Conic"],
PARAMETER["False_Easting",6561666.66666667],PARAMETER["False_Northing",1640416.66666667],
PARAMETER["Central_Meridian",-120.5],PARAMETER["Standard_Parallel_1",37.0666666666667],
PARAMETER["Standard_Parallel_2",38.4333333333333],PARAMETER["Latitude_Of_Origin",36.5],UNIT["Foot_US",0.3048006096012192]]

SELECT srid, srtext, proj4text 
FROM spatial_ref_sys 
    WHERE  srtext ILIKE '%California%' 
        AND proj4text ILIKE '%ft%' AND srtext ILIKE '%zone 3%';

shp2pgsql: Load data

Loads into Geography format

shp2pgsql -D -I -s 4326 -G  ne_10m_populated_places data.world_cities | psql

Loads into geometry

shp2pgsql -D -I -s 4269  tl_2014_us_state data.us | psql
shp2pgsql -D -I -t 2D -s 4269  BART_13 data.bart | psql -d workshop_1
shp2pgsql -D -I -t 2D -s 4269  BART_Sta_13 data.bart_stations | psql -d workshop_1
shp2pgsql -D -I -s 2227 building_footprint data.sfo_buildings | psql
shp2pgsql -D -I -s 4326 -t 2D london_tube_lines data.london_tube_lines | psql

Using ogr2ogr to load data

If you started with the raw KML file, you can do this from command line. On windows get rid of \ and put all on one line

ogr2ogr -f "PostgreSQL"  \
PG:"host=localhost user=postgres port=5432 dbname=workshop_1 password=whatever"  \
"London tube Lines.kml"  -nln data.london_tube_lines -lco GEOMETRY_NAME=geom

3D Examples

Build 3D buildings by Extrusion

DROP TABLE IF EXISTS data.sfo_3dbuildings; 
-- 3477 rows, 45763 ms
SELECT gid, objname, (g).path[1] As poly_id, 
    ST_Translate(
        ST_Extrude(
            ST_MakePolygon( ST_ExteriorRing((g).geom) ) ,0,0,z_max - z_min)
    ,0,0,z_min)::geometry(POLYHEDRALSURFACEZ,2227) As geom
INTO data.sfo_3dbuildings
FROM (SELECT gid, objname, ST_Dump(geom) AS g, z_max, z_min
FROM data.sfo_buildings 
WHERE ST_Expand(geom,5000 ) 
    && ST_Transform(ST_SetSRID(ST_Point(-122.404143, 37.79322),4326),2227)
) As f;

Add 3D Index

CREATE INDEX idx_sfo_3dbuilds_geom
   ON data.sfo_3dbuildings USING gist(geom gist_geometry_ops_nd);

X3D query, randomly color buildings

Use 3D bounding box &&& operator and form a 3D box filter

SELECT string_agg('<Shape><Appearance>
  <Material ambientintensity="0.200" 
    containerfield="material" shininess="0" diffusecolor="0.5 ' 
        || random()::numeric(6,4)::text || ' ' 
        || random()::numeric(6,4)::text || '"/>
  </Appearance>' 
    || ST_AsX3D(geom)  || '</Shape>', '')
FROM data.sfo_3dbuildings
WHERE geom &&& ST_Expand(
    ST_Force3D(
        ST_Transform(
            ST_SetSRID(ST_Point(-122.404143, 37.79322),4326),2227)), 1000);

View in X3Dom

X3D query, Using texture

Use 3D bounding box &&& operator and form a 3D box filter

SELECT string_agg('<Shape><Appearance>
  <ImageTexture url=''"images/bt'  || 
    CASE WHEN 
        (ST_XMax(geom) - ST_XMin(geom)) > (ST_Zmax(geom) - ST_ZMin(geom)) 
        THEN ceiling(random()*4)::integer::text 
    WHEN ST_ZMax(geom) < 90 THEN ceiling(3 + random()*3)::integer::text
ELSE ceiling(5 + random()*8)::integer::text END
  || '.jpg"'' />
  </Appearance>' 
    || ST_AsX3D(geom)  || '</Shape>', '')
FROM data.sfo_3dbuildings
WHERE geom &&& ST_Expand(
    ST_Force3D(
        ST_Transform(
            ST_SetSRID(ST_Point(-122.404143, 37.79322),4326),2227)), 2000);

X3Dom with texture

X3D query, Draw drone path

Launch a drone and graph out its trajectory

WITH drone_path AS (
    SELECT  ST_GeomFromText('LINESTRING Z (6012500 2121900 0, 
    6011608 2116479 100,
    6011915 2117592 200,6012569 2116636 150,
    6013788 2115817 90, 6013790 2115830 0)',2227)  As geom)
SELECT '<Shape>' 
    || ST_AsX3D(geom)  || 
    '<Appearance>
  <Material  containerfield="material"  emissiveColor="0.7 0.2 0.2" />
  </Appearance></Shape>'
FROM drone_path;

X3D query, Highlight buildings in red on drone path

Use ST_3DIntersects buildings intersect path

set postgis.backend=geos;
WITH drone_path AS (SELECT 
 ST_GeomFromText('LINESTRING Z (6012500 2121900 0, 
    6011608 2116479 100,
    6011915 2117592 200,6012569 2116636 150,
    6013788 2115817 90, 6013790 2115830 0)',2227) As geom)
SELECT string_agg('<Shape><Appearance>
 <Material ambientintensity="0.200" containerfield="material" 
    transparency="0.1" shininess="0" diffusecolor="1 1 0"/>
  </Appearance>' 
    || ST_AsX3D(b.geom)  || '</Shape>', '')
FROM data.sfo_3dbuildings As b 
    INNER JOIN drone_path As d ON ST_3DIntersects(b.geom, d.geom);

Drone trajectory buildings in red