View source for Loading data into PostGIS

This page is now getting a little more attention and organization.

Refer to "PostGIS in Action" (Obe and Hsu) Chapter 7

I have written a script that handles loading shapefiles, geodatabases, CSV tables, DBF tables, and USDA SSURGO data into PostGIS.

== Loading shapefiles ==

(Here is my quick ref of [[EPSG projection codes]])

=== First method: Loading shapefiles with shp2pgsql ===

Use the shp2pgsql command to create and load the table and then manually define the projection. Example:

 shp2pgsql roads.shp clatsop.roads | psql -U gis_owner -d gis_data 

The shp2pgsql command appears to ignore the PRJ file that defines projection. The data is loaded as projected
but the SRID column will be NULL, so I manually enter a command to define the projection like this:

 psql -U gis_owner -d gis_data
 SELECT UpdateGeometrySRID('clatsop.roads', 'the_geom', 2913);

There is also a plugin called shp2pgsql-gui included with the PostGIS source code for pgAdmin III that lets you load shapefiles directly from the pgAdmin tool. This looks elegant but I have not tested it yet.

=== Second method: Loading almost any data with ogr2ogr ===

I have information about using 
''ogr2ogr'' to load GPX data into PostGIS on the [[GPX]] page.

Using '''ogr2ogr'' is convenient because you can load the data and reproject it in one step. Note source and destination are reversed in this command.
Also very convenient because you can load any data that ogr2ogr knows about, do '''ogr2ogr --formats''' for the list,
not just shapefiles.

 ogr2ogr -skipfailures -t_srs EPSG:2913 -f PostgreSQL PG:"host=localhost user=gis_owner dbname=gis_data" \
 roads.shp roads \
 -lco GEOMETRY_NAME=the_geom -lco SCHEMA=clatsop -lco OVERWRITE=YES \
 -nlt MULTILINESTRING

I chose to force the geometry type to "MULTI" to avoid conversion failures; the alternative would be to convert
the shapefile to simple line features first.

Do note the order is "destination source" not "source destination"

For [[Mapping Vietnam]] I want to load data from GIS Data Depot into PostGIS. At the moment I only want to do this so I have some data for testing. Later on I might even want these layers in a map.

The downloaded data are in ESRI "E00" format and have been compressed.

 # Uncompress all downloaded files
 for file in *.E00.GZ; do gunzip $file; done
 
 # put the data into "Web Mercator" projection
 PRJ="EPSG:900913"
 DB=vietnam

 # The original files were coverages, so I have to know something about
 # what I want to transfer. For example for this file I want only PAL
 '''ogrinfo PPPOLY.E00'''
 1: ARC (Line String)
 2: CNT (Point)
 3: LAB (Point)
 ''4: PAL (Polygon)'' This is the one I want.
 
 ogr2ogr -t_srs $PRJ -f PostgreSQL \
 PG:"user=postgres dbname=$DB" \
 PPPOLY.E00 PAL \
 -lco GEOMETRY_NAME=the_geom -nln political_boundaries

== Loading broken shapefiles ==

You could simply drop all the broken bits and pieces,
or you could try to repair the uncooperative rows.

I am beginning to think that ''all'' shapefiles are "broken" and that this fact might be why ArcGIS models crash randomly.

'''I jumped through many hoops (see below for now)
trying to get shapefiles to load with ogr2ogr.
Then I found that shp2pgsql just worked.
'''

 shp2pgsql -c -s 32026 $SRCDIR/centerline.shp > centerline.sql
 psql -U postgres -f centerline.sql solarmap

=== Dropping broken rows ===

First off, if you have access to ArcGIS, run your shapefiles through the
"Repair Geometry" tool to repair or remove broken features.

Secondly, apparently it is common for shapefiles to have columns that exceed the constraints defined in their DBF fork. Use PRECISION=NO to force them to load. 

Thirdly, use 'skipfailures' to drop features that are wrong, for example in multilinestring format in a linestring shapefile.

 ogr2ogr -skipfailures -t_srs $PRJ $DST $SRC -lco GEOMETRY_NAME=the_geom -lco PRECISION=NO

=== Repairing ===

BUT WAIT! I really don't want to simply drop the benighted features. I want to fix them. And by the way, I don't happen to have ArcGIS right now. 

First generate a SQL file.

 '''ogr2ogr -t_srs $PRJ -f PGDump $SRC.sql $SRC -lco GEOMETRY_NAME=the_geom -lco PRECISION=NO'''

Next edit the file to force the AddGeometryColumn function to create a GEOMETRY type column instead of a LINESTRING column. This will allow any combination of LINESTRING and MULTILINESTRING objects to be loaded.

Finally, load the SQL data using psql

 '''psql -U postgres -d corvallis -f $SRC.sql'''

This worked; I now have a centerlines file with 5218 rows in it. See?

 corvallis=# '''select COUNT(ogc_fid) from centerline;'''
 count 
 -------
  5218

Here are 2 queries that list the MULTISTRING rows

 SELECT ogc_fid,fullname FROM centerline WHERE GeometryType(geom_nad27)='MULTILINESTRING';

 SELECT ogc_fid,fullname FROM centerline WHERE ST_NumGeometries(geom_nad27)>1;

Break into pieces (works for single part objects too)

 SELECT ogc_fid,ST_Dump(geom_nad27) FROM centerline WHERE ST_NumGeometries(geom_nad27)>1;

=== Comparison ===

Here are commands I used to load the data for the [[ArcGIS For Server]] comparison article.

First I need a database.

Then I need data, loaded.

=== Check your work ===

The ogrinfo command can peek into your PostGIS database just as it can look into any spatial file it knows how to read. The -so option gives summary output instead of dumping every row to the screen.

 ogrinfo -so PG:"user=postgres dbname=vietnam" political_boundaries

== Loading data from ESRI File Geodatabases ==

You can load data without the ESRI API now just using the OpenFGDB driver built into GDAL.

If the files are in 10.x format then you can use GDAL. 
If you have old files (< v10) you can update them first using ArcGIS ArcCatalog -- use "Properties" and click "Update" button.
Before you start make sure the FGDB is not locked (for example, open in a map) and that you have write permission. If you don't do this then the tool will grind away for perhaps 20 minutes before telling you something is wrong.

I created a database to hold the new tables.

 createdb -U postgres or_taxlots -T template_postgis

Next I listed of all tables that live in this FGDB. It looks like this file has taxlots polygons and attributes for just about the whole state.
If I wanted to import the entire file I'd use a script.

 ogrinfo Taxlots.gdb

The Metro data is in: NAD_1983_HARN_StatePlane_Oregon_North_FIPS_3601 which is EPSG:2913 see [[EPSG projection codes]]

This page helped me a lot: http://www.postgresonline.com/article_pfriendly/31.html

Just working on metro at the moment so I am not going to create the script right now to load everything.

# This generates a file containing SQL, as a test but we could load it too using psql command.
# The ogr2ogr command will make the table name lowercase for us -- taxlots_metro instead of the illegal taxlots_METRO.

 ogr2ogr \
 -s_srs "EPSG:2913" \
 -t_srs "EPSG:2994" -f PostgreSQL PG:"host=dart user=gis_owner dbname=gis_data" \
 Taxlots.gdb taxlots_METRO  \
 -lco GEOMETRY_NAME=geom_ogic -lco OVERWRITE=yes -lco SCHEMA=or_state -lco GEOMETRY_NAME=geom

This particular set has lots of attribute data and I'd like to work on a script to clean up the data as I import it.

 ogr2ogr -s_srs "EPSG:2992" -t_srs "EPSG:2994" $DST orcnty24.shp orcnty24 -lco GEOMETRY_NAME=geom
 ogr2ogr -s_srs "EPSG:2992" -t_srs "EPSG:2994" $DST or_state_boundary.shp or_state_boundary -lco GEOMETRY_NAME=the_geom

This worked for me uploading GNIS data from a FGDB to PostGIS:

 ogr2ogr -f PostgreSQL PG:'dbname=gis_data' -nln or_state.gnis GNIS_OR.gdb GNIS_OR_Buffer

I loaded this as postgres then changed ownership to gis_data. It would not let me load it as gis_owner. (?)

== Loading GNIS data from USGS BGN ==

Download data from http://geonames.usgs.gov/domestic/download_data.htm

For the whole country:
 wget http://geonames.usgs.gov/docs/stategaz/NationalFile_20130602.zip

Convert pipes to tabs
 zcat NationalFile_20130602.zip | sed  's/|/\t/g' > gnis.csv

Create target schema
 psql -c "CREATE SCHEMA usgs;" -U gis_owner gis_data

Create a gnis.vrt file containing this

 cat > gnis.vrt
 <OGRVRTDataSource>
    <OGRVRTLayer name="gnis">
        <SrcDataSource>gnis.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="PRIM_LONG_DEC" y="PRIM_LAT_DEC"/>
    </OGRVRTLayer>
 </OGRVRTDataSource>

Load!
 ogr2ogr -t_srs EPSG:4326 -s_srs EPSG:4326 -f PostgreSQL PG:"user=gis_owner dbname=gisdata active_schema=usgs" gnis.vrt -nln national_gnis gnis 

Get rid of the now redundant csv file
 rm gnis.csv

== Loading census data from TXT files ==

This is almost the same as loading the BGN GNIS data.
Here is a script to unzip and load a series of Gaz* files.

 psql -c "CREATE SCHEMA us_census;" -U gis_owner gisdata
 export PGCLIENTENCODING="latin1"

 # Import Gazetteers
 for zip in G*.zip; do 
  base=`basename $zip .zip`
  echo $base

  zcat $zip > $base.csv
  cat <<EOF > $base.vrt
 <OGRVRTDataSource>
   <OGRVRTLayer name="$base">
       <SrcDataSource>$base.csv</SrcDataSource>
       <GeometryType>wkbPoint</GeometryType>
       <LayerSRS>WGS84</LayerSRS>
       <GeometryField encoding="PointFromColumns" x="INTPTLONG" y="INTPTLAT"/>
   </OGRVRTLayer>
 </OGRVRTDataSource>
 EOF
  #ogrinfo -so $base.vrt $base
  psql -c "DROP TABLE us_census.$base" -U gis_owner gisdata
  ogr2ogr -a_srs EPSG:4326 -f PostgreSQL PG:"user=gis_owner dbname=gisdata active_schema=us_census" $base.vrt -nln $base $base -lco SEPARATOR=TAB -lco ENCODING=LATIN1

 rm $base.csv $base.vrt
 done

== Loading data from KMZ files ==

I exported some data from ArcMap into KMZ files (using "Map to KML" which worked, "Layer to KML" did not.)

Now I want to load them into PostGIS.

Next I want to make the data available via Geoserver and define styles there, too.

GDAL can't handle compressed (KMZ) files so first I have to unpack them.

 for i in /t/Brian/Android/*kmz; do name=`basename $i .kmz`; mkdir $name; cd $name; unzip $i; cd ..; done 

Now I can use ogrinfo to examine each doc.kml file. 

 for i in *; do ogrinfo -ro $i/doc.kml; done 

 INFO: Open of `corvallis_buildings/doc.kml' 
      using driver `KML' successful.
 1: Buildings (Corvallis) (Multi Polygon)
 INFO: Open of `corvallis_fire_symbols/doc.kml'
      using driver `KML' successful.
 1: Corvallis Fire Symbols (Point)
 INFO: Open of `corvallis_parks/doc.kml'
      using driver `KML' successful.
 1: Corvallis parks (Multi Polygon)
 INFO: Open of `corvallis_trails/doc.kml'
      using driver `KML' successful.
 1: Trails (Corvallis parks) (Multi Line String)
 INFO: Open of `oregon_emergency_response/doc.kml'
      using driver `KML' successful.
 1: Emergency Response (statewide) (Point)
 INFO: Open of `or_linnbentonmarion_roads/doc.kml'
      using driver `KML' successful.
 1: Main Routes (Multi Line String)

Load each data set, thusly

 cd corvallis_buildings
 ogr2ogr -f PostgreSQL PG:"user=postgres dbname=hydrant_survey" doc.kml  'Buildings (Corvallis)'
 cd ..

I used pgadminIII to rename each resultant table to remove any stupid characters (eg spaces) from its name.

Now that I have all data loaded into PostGIS, I preview it with QGIS to make sure it still looks usable.

Next I go into Geoserver and create stores and layers for each new dataset.

== Loading data from OpenStreetMap ==

Moved to its own page [[OpenStreetMap]]

== Loading VMAP data ==

I got help here: http://tim.linfiniti.com/post/68195811/using-vmap-data-with-gdal-ogr

Download and unpack in a directory tree that has '''no upper case letters'''. Be careful about '''trailing slashes''', This driver is touchy. I needed Vietnam so I got the south Asia and Australia dataset.

To list layers, 
 ogrinfo -ro -so '''gltp:/vrf'''/home/bwilson/vmap_raw_vpf/v0sas_5/vmaplv0/sasaus

To load layers into PostGIS,
 ogr2ogr -a_srs EPSG:4326 -f PostgreSQL PG:dbname=${DB} -update -append -nln ${LAYERBASE} gltp:/vrf/${DATAPATH} �${LAYER}'

=== Check your work ===

 ogrinfo -so PG:"user=postgres dbname=vmap" '''layername'''

== Reprojecting ==

Now that I have some data loaded I need to put it into the projection I need.

== In the ESRI world ==

* I want to be able to move data from ESRI proprietary formats into PostGIS.
* I want to be able to hand a set of ESRI compatible tools to my GIS Analyst co-workers so they can get the data into my fabled PostGIS data warehouse without my help. (Using ESRI "toolboxes" and ESRI Model Builder for example)

Note although I might be able to use an ArcSDE license (aka ArcGIS Server) to access PostGIS but I regard that as so much extra work it's not worth it. Also it requires an ArcInfo license on the desktop and [http://www.bartleby.com/129/ I would prefer not to.]

You only need an ArcGIS license for this to access the ESRI formats. If your data are in shapefiles you don't need this... scroll up.

In fact if your data are in ESRI File Geodatabases that are new (version 10) you can use GDAL 1.9.0 to access them, I can do this with GDAL on Linux, even. But most of my FGDBs are older. 

In fact you only need the license for a few things these days. 

=== Iterate a feature class ===

Here is a simple script (9.3) to iterate over a feature class and read its geometry.
The nice thing about using the ESRI code for this is that it does not matter what the data source is, it's just a 'feature class' that can be stored in a shapefile or a personal geodatabase or ArcSDE... etc...

<pre>
#!/usr/bin/env python

import arcgisscripting 
gp = arcgisscripting.create(9.3)

# Hard coded data source, for simplicity.
workspace = 'D:/AGIProducts/IncidentView_Data/Data/TEMPORARY_WORKSPACE/WA_King.gdb'
featureclass = 'test_Addresses_points'

gp.workspace = workspace
desc = gp.describe(featureclass) # Get some metadata
shapefieldname = desc.ShapeFieldName
print "ShapeType %s" % desc.ShapeType
print "Indexed =", desc.HasSpatialIndex
print "Shape field", shapefieldname

fields = gp.ListFields(featureclass, '*') # Get a list of the attributes

rows = gp.SearchCursor(featureclass)
row = rows.Next()
i = 0
while row:
    feature = row.GetValue(shapefieldname)
    # If this is a point, we can just grab its shape
    point = feature.GetPart()
    print i, point.x, point.y
    row = rows.Next()
    i += 1

del gp # Release the geoprocessing object
# That's all!
</pre>

==== Write the feature to PostGIS ====

So now I want to put the points into PostGIS... how to do that?

I probably could use the GDAL python bindings.
I can probably figure out how to write the geometry directly to PostgreSQL

== Loading Raster data ==

=== Hillshade ===

I have an existing Oregon hillshade in a TIFF file.
I don't want to do any analysis with it so this is probably not a good idea.
I will simply make it available for now via NFS file share so I can view it in QGIS.

Later I will use it as a background in a tiled map.

== References ==

Paulo Corti's notes on
[http://www.paolocorti.net/2006/07/27/migrating-shapefiles-to-postgis/ Migrating shapefiles to PostGIS] from "Thinking in GIS"

[[Category: GIS|PostGIS]]
[[Category: PostgreSQL]]