Paragon Corpoation PostGIS Spatial Database Engine OSGeo.org The Open Source Geospatial Foundation UMN Mapserver Boston Geographic Information Systems   
FOSS4G International 2017, August 14th-18th 2017
   PostGreSQL Object Relational Database Management System
Home   About Boston GIS   Consulting Services  Boston GIS Blog  Postgres OnLine Journal  Planet PostGIS  PostGIS Funding

Purpose of BostonGIS

BostonGIS is a testbed for GIS and Web Mapping solutions utilizing open source, freely available and/or open gis technologies. We will be using mostly Boston, Massachusetts data to provide mapping and spatial database examples.

If you have some thoughts or comments on what you would like to see covered on this site, drop us a line on our Feed Back page.


GIS Tutorials on Opensource and OpenGIS technologies Tutorials
GIS Article comments Article and Tutorial Comments
Boston GIS BLog Rss FeedBoston GIS blog

PDF HTML All BostonGIS tutorials packaged together in an E-Book.


Boston GIS Store

Loading


Tutorial and Tip Sites
Desktop GIS
External Data
GIS Events and Groups
GIS SDKs and Frameworks
External Resources
Glossary
GIS Blogs Around Boston
External GIS Blogs
External Papers Articles
GIS Quick Guides and References
OpenStreetMap and OpenLayers Tutorials
PostGIS, pgRouting, and PostgreSQL Tutorials
Part 1: Getting Started With PostGIS: An almost Idiot's Guide (PostGIS 2.2) more ...
pgRouting: Loading OpenStreetMap with Osm2Po and route querying more ...
Part 1: Getting Started With PostGIS: An almost Idiot's Guide (PostGIS 2.0) more ...
OSCON 2009: Tips and Tricks for Writing PostGIS Spatial Queries more ...
PGCon2009: PostGIS 1.4, PostgreSQL 8.4 Spatial Analysis Queries, Building Geometries, Open Jump more ...
PLR Part 3: PL/R and Geospatial Data Abstraction Library (GDAL) RGDAL more ...
PostGIS Nearest Neighbor: A Generic Solution - Much Faster than Previous Solution more ...
Solving the Nearest Neighbor Problem in PostGIS more ...
PLR Part 2: PL/R and PostGIS more ...
PLR Part 1: Up and Running with PL/R (PLR) in PostgreSQL: An almost Idiot's Guide more ...
Part 2 - PostGIS and SharpMap in ASP.NET 2.0 using VB.NET: Displaying the Maps more ... download
Part 1 - PostGIS and SharpMap in ASP.NET 2.0 using VB.NET: Compiling SharpMap with PostGIS more ...
Part 3: PostGIS Loading Data from Non-Spatial Sources

Printer Friendly

Often you will receive data in a non-spatial form such as comma delimited data with latitude and longitude fields. To take full advantage of PostGIS spatial abilities, you will want to create geometry fields in your new table and update that field using the longitude latitude fields you have available.

General Note: All the command statements that follow should be run from the PgAdminIII Tools - Query Tool or any other PostGreSQL Administrative tool you have available. If you are a command line freak - you can use the psql command line tool packaged with PostGreSQL.

Getting the data

For this exercise, we will use US zip code tabulation areas instead of just Boston data. The techniques here will apply to any data you get actually.

First step is to download the data from US Census. http://www.census.gov/geo/www/gazetteer/places2k.html

Importing the Data into PostGreSQL

PostgreSQL comes with a COPY function that allows you to import data from a delimited text file. Since the ZCTAs data is provided in fixed-width format, we can't import it easily without first converting it to a delimited such as the default tab-delimited format that COPY works with. Similarly for data in other formats such as DBF, you'll either want to convert it to delimited using tools such as excel, use a third party tool that will import from one format to another, or one of my favorite tools Microsoft Access that allows you to link any tables or do a straight import and export to any ODBC compliant database such as PostgreSQL.

For fixed width data such as ZCTA, we demonstrate a more generic way of importing fixed with data into PostgreSQL that doesn't require any additional tools beyond what is packaged with PostgreSQL.

Create the table to import to

First you will need to create the table in Postgres. You want to make sure the order of the fields is in the same order as the data you are importing.


CREATE TABLE zctas
(
  state char(2),
  zcta char(5),
  junk varchar(100),
  population_tot int8,
  housing_tot int8,
  water_area_meter float8,
  land_area_meter float8,
  water_area_mile float8,
  land_area_mile float8,
  latitude float8,
  longitude float8
) 
WITHOUT OIDS;

Convert from Fixed-width to Tab-Delimited

For this part of the exercise, I'm going to use Microsoft Excel because it has a nice wizard for dealing with fixed-width and a lot of windows users have it already. If you open the zcta file in Excel, it should launch the Text Import Wizard. MS Access has a similarly nice wizard and can deal with files larger than excels 65000 some odd limitation. Note there are trillions of ways to do this step so I'm not going to bother going over the other ways. For non-MS Office users other office suites such as Open-Office probably have similar functionality.

  1. Open the file in Excel.
  2. Import Text Wizard should launch automatically and have Fixed-Width as an option
  3. Look at the ZCTA table layout spec http://www.census.gov/geo/www/gazetteer/places2k.html#zcta and set your breakouts the same as specified. For the above I broke out the Name field further into first 5 for zcta and the rest for a junk field.
  4. Next File->Save As ->Text (Tab delimited)(*.txt) -give it name of zcta5.tab
  5. Copy the file to somewhere on your PostGreSQL server.
  6. The COPY command

    Now copy the data into the table using the COPY command. Note the Copy command works using the PostGreSQL service so the file location must be specified relative to the Server.

    
         COPY zctas FROM 'C:/Downloads/GISData/zcta5.tab';
    
    

    Creating and Populating the Geometry Field

    Create the Geometry Field

    To create the Geometry field, use the AddGeometryColumn opengis function. This will add a geometry field to the specified table as well as adding a record to the geometry_columns meta table and creating useful constraints on the new field. A summary of the function can be found here http://postgis.refractions.net/docs/ch06.html#id2526109.

    SELECT AddGeometryColumn( 'public', 'zctas', 'thepoint_lonlat', 4269, 'POINT', 2 );

    The above code will create a geometry column named thepoint_longlat in the table zctas that validates to make sure the inputs are 2-dimensional points in SRID 4269 (NAD83 longlat).

    Populate the Geometry Field using the Longitude and Latitude fields

    
    UPDATE zctas
            SET thepoint_lonlat = ST_SetSRID(ST_Point( longitude, latitude),4269) )
    
    

    The above code will generate a PostGIS point geometry object of spatial reference SRID 4269.

    There are a couple of things I would like to point out that may not be apparently clear to people not familiar with PostGreSQL or PostGis

    • You can't just put any arbitrary SRID in there and expect the system to magically transform to that. The SRID you specify has to be the reference system that your data is in.

    Transforming to Another spatial reference system

    The above is great if you want your geometry in longlat spatial reference system. In many cases, longlat is not terribly useful. For example if you want to do distance queries with your data, you don't want your distance returned back in longlat. You want it in a metric that you normally measure things in.

    In the code below, we will create a new geometry field that holds points in the WGS 84 North Meter reference system and then updates that field accordingly.

    
    SELECT AddGeometryColumn( 'public', 'zctas', 'thepoint_meter', 32661, 'POINT', 2 );
    
    UPDATE zctas SET thepoint_meter = transform(setsrid(makepoint(longitude, latitude),4269), 32661);
    
    

    Please note in earlier versions of this article I had suggested doing the above with

    
    UPDATE zctas
    SET thepoint_meter = transform(PointFromText('POINT(' || longitude || ' ' || latitude || ')',4269),32661) ;
    
    

    On further experimentation, it turns out that PointFromText is perhaps the slowest way of doing this in Postgis. Using a combination of ST_Setsrid and ST_Point is on the magnitude of 7 to 10 times faster at least for versions of PostGIS 1.2 and above. ST_GeomFromText comes in second (replace ST_PointFromText with ST_GeomFromText) at about 3 to 1.5 times slower than ST_SetSRID ST_Point. See about PointFromText on why PointFromText is probably slower. In ST_GeomFromText, there appears to be some caching effects so on first run with similar datasets ST_GeomFromText starts out about 3-4 times slower than the ST_makepoint (ST_Point) way and then catches up to 1.5 times slower. This I tested on a dataset of about 150,000 records and all took - 26 secs for ST_PointFromText fairly consistently, 10.7 secs for first run of GeomFromText then 4.1 secs for each consecutive run, 3.5 secs fairly consistently for setsrid,makepoint on a dual Xeon 2.8 Ghz, Windows 2003 32-bit with 2 gig RAM).

    Index your spatial fields

    One of the number one reasons for poor query performance is lack of attention to indexes. Putting in an index can make as much as a 100 fold difference in query speed depending on how many records you have in the table. For large updates and imports, you should put your indexes in after the load, because while indexes help query speed, updates against indexed fields can be very slow because they need to create index records for the updated/inserted data. In the below, we will be putting in GIST indexes against our spatial fields.

    
    CREATE INDEX idx_zctas_thepoint_lonlat ON zctas
      USING GIST (thepoint_lonlat);
    
    CREATE INDEX idx_zctas_thepoint_meter ON zctas
      USING GIST (thepoint_meter);
    
    ALTER TABLE zctas ALTER COLUMN thepoint_meter SET NOT NULL;
    CLUSTER idx_zctas_thepoint_meter ON zctas;
    
    VACUUM ANALYZE zctas;
    
    

    In the above after we create the indexes, we put in a constraint to not allow nulls in the thepoint_meter field. The not null constraint is required for clustering since as of now, clustering is not allowed on gist indexes that have null values. Next we cluster on this index. Clustering basically physically reorders the table in the order of the index. In general spatial queries are much slower than attribute based queries, so if you do a fair amount of spatial queries especially bounding box queries, you get a significant gain.

    Please note: (for those familiar with Cluster in SQL Server - Cluster in PostgreSQL pretty much means the same thing, except in PostgreSQL - at least for versions below 8.3 (not sure about 8.3 and above), if you cluster and then add new data or update data, the system does not ensure that your new or updated data is ordered according to the cluster as it does in SQL Server. The upside of this is that there is no additional penalty in production mode with clustering as you have in SQL Server. The downside is you must create a job of some sort to maintain the cluster for tables that are constantly updated.

    In the above we vacuum analyze the table to insure that index statistics are updated for our table.





    Post Comments About Part 3: PostGIS Loading Data from Non-Spatial Sources
Part 2: Introduction to Spatial Queries and SFSQL with PostGIS more ...
Miscellaneous Tutorials/Cheatsheets/Examples
SpatiaLite Tutorials
Boston External Map Examples
SQL Server 2008 Tutorials
UMN Mapserver Tutorials
General Commentary
 CommenterComment
12/23/2009 12:45:29 PMReginaBalu Nadig,

You are using ST_MakePoint wrong. The third argument is not the SRID, but the Z coordinate if you have one (there are 2 ST_MakePoint variants one that takes X,Y and one that Takes X,Y,Z).

You want to do

ST_SetSRID(ST_MakePoint(longitude,latitude), yoursridhere)
12/18/2009 12:27:16 PMBalu NadigThnaks for your very helpful post...

I need some help... I have created a postgis database and imported borough and districts data. I am trying do a google map of organizations (have their lat/lng in the database) within the borough/districts... I created boros with srid 2831... Here is my problem...

1)
select boroname from boros b
where ST_Contains(setsrid(b.the_geom, 2831) , st_makepoint(40.7195767, -73.9857889, 2831))

ERROR: Operation on two geometries with different SRIDs

2)
select boroname from boros b
where ST_Contains((b.the_geom) , st_makepoint(40.7195767, -73.9857889, 2831))

ERROR: Operation on two geometries with different SRIDs

3)
select boroname from boros b
where ST_Contains(setsrid(b.the_geom, -1) , st_makepoint(40.780704, -73.972342,-1))

No error, but no results either :(

Thanks for your help :)
11/11/2009 10:20:37 AMReginaIt is expressed in degrees. USe ST_Distance_Sphere to get meters (if you only have points), and use an ST_Transform to something like 2163 to get approximate meters. for versions prior to 1.5, ST_Distance_Sphere only works with points, but from postGIS 1.5 - ST_Distance_Sphere/Spheroid works with all geometries except for curved geometries.
11/4/2009 3:59:54 AMChrisHi,
I have loaded the data, and set the SRID.

Can you explain me results (3.57...) following query:

select st_distance(ny.thepoint_lonlat, wa.thepoint_lonlat)
from zctas ny, zctas wa
where ny.junk = 'New York city'
and (wa.junk = 'Washington city' and wa.state = 'DC');

It is expressed in ... ?
10/22/2009 9:20:46 PMLeoThanks. Fixed
10/22/2009 3:24:25 PMedwardIthink ST_Porint should be ST_Point?
9/25/2009 2:56:00 AMReginaThomas,

Both 32661 and 4326 use the WGS 84 datum and ellipsoid but one is long lat (4326) and one is a meter based one. The WGS 84 simply defines the ellipse and datum.
If you look closely -- We are transforming from 4269 (which is pretty much 4326 WGS 84 (its NAD 80) long lat) to 32661 a meter based projection more suitable for measurment. If you query the PostGIS spatial_ref_sys table

SELECT srid, srtext, proj4text
FROM spatial_ref_sys where srid = 32661;

you'll see the detailed info about the spatial ref.

I think spatial reference systems are the most confusing thing in GIS and tough to explain.
9/21/2009 11:47:29 AMthomasI just wanted to know what the difference between SRID '32661' and '4326'? I thought SRID 4326 is the identifier for the WGS84 reference system. But in this example you use SRID '32661' for WGS84.
9/1/2009 9:46:49 AMNadHi,

I am trying to build a trigger that implements automatic population of geometry fields in my table. Following the example with table zctas, i wrote the following;
--
CREATE OR REPLACE Function zctaspoint_conv() RETURNS TRIGGER AS
$$

BEGIN
UPDATE zctas SET thepoint_lonlat = PointFromText('POINT(' || longitude || ' ' || latitude || ')',4269);

END;
$$ LANGUAGE 'plpgsql';

CREATE TRIGGER zctaspoint_conv
AFTER INSERT OR UPDATE
ON zctas FOR EACH ROW
EXECUTE PROCEDURE zctaspoint_conv();
--

however, the trigger fails when i do inserting or updating in zctas table. What did I miss?
Thanks for any feedback.

nad
Locations of visitors to BostonGIS
Boston GIS      Copyright 2017      Paragon Corporation