PostGIS Cookbook - Second Edition

We will import the firenews.csv file that stores a series of web news collected from various RSS feeds related to forest fires in Europe in the context of the European Forest Fire Information System (EFFIS), available at http://effis.jrc.ec.europa.eu/.

For each news feed, there are attributes such as place name, size of the fire in hectares, URL, and so on. Most importantly, there are the x and y fields that give the position of the geolocalized news in decimal degrees (in the WGS 84 spatial reference system, SRID = 4326).

For Windows machines, it is necessary to install OSGeo4W, a set of open source geographical libraries that will allow the manipulation of the datasets. The link is: https://trac.osgeo.org/osgeo4w/

In addition, include the OSGeo4W and the Postgres binary folders in the Path environment variable to be able to execute the commands from any location in your PC.

The steps you need to follow to complete this recipe are as shown:

      $ cd ~/postgis_cookbook/data/chp01/
      $ head -n 5 firenews.csv

The output of the preceding command is as shown:

      $ psql -U me -d postgis_cookbook 
      postgis_cookbook=> CREATE EXTENSION postgis; 
      postgis_cookbook=> CREATE SCHEMA chp01; 
      postgis_cookbook=> CREATE TABLE chp01.firenews 
      ( 
        x float8, 
        y float8, 
        place varchar(100), 
        size float8, 
        update date, 
        startdate date, 
        enddate date, 
        title varchar(255), 
        url varchar(255), 
        the_geom geometry(POINT, 4326) 
      );

      postgis_cookbook=> COPY chp01.firenews (
        x, y, place, size, update, startdate, 
        enddate, title, url
      ) FROM '/tmp/firenews.csv' WITH CSV HEADER;

      postgis_cookbook=> SELECT COUNT(*) FROM chp01.firenews;

The output of the preceding command is as follows:

     postgis_cookbook=# SELECT f_table_name, 
     f_geometry_column, coord_dimension, srid, type 
     FROM geometry_columns where f_table_name = 'firenews';

The output of the preceding command is as follows:

      postgis_cookbook=> UPDATE chp01.firenews 
      SET the_geom = ST_SetSRID(ST_MakePoint(x,y), 4326); 
      postgis_cookbook=> UPDATE chp01.firenews 
      SET the_geom = ST_PointFromText('POINT(' || x || ' ' || y || ')',
                                      4326);

      postgis_cookbook=# SELECT place, ST_AsText(the_geom) AS wkt_geom 
      FROM chp01.firenews ORDER BY place LIMIT 5;

The output of the preceding comment is as follows:

      postgis_cookbook=> CREATE INDEX idx_firenews_geom
      ON chp01.firenews USING GIST (the_geom);

This recipe showed you how to load nonspatial tabular data (in CSV format) in PostGIS using the COPY PostgreSQL command.

After creating the table and copying the CSV file rows to the PostgreSQL table, you updated the geometric column using one of the geometry constructor functions that PostGIS provides (ST_MakePoint and ST_PointFromText for bi-dimensional points).

These geometry constructors (in this case, ST_MakePoint and ST_PointFromText) must always provide the spatial reference system identifier (SRID) together with the point coordinates to define the point geometry.

Each geometric field added in any table in the database is tracked with a record in the geometry_columns PostGIS metadata view. In the previous PostGIS version (< 2.0), the geometry_fields view was a table and needed to be manually updated, possibly with the convenient AddGeometryColumn function.

For the same reason, to maintain the updated geometry_columns view when dropping a geometry column or removing a spatial table in the previous PostGIS versions, there were the DropGeometryColumn and DropGeometryTable functions. With PostGIS 2.0 and newer, you don't need to use these functions any more, but you can safely remove the column or the table with the standard ALTER TABLE, DROP COLUMN, and DROP TABLE SQL commands.

In the last step of the recipe, you have created a spatial index on the table to improve performance. Please be aware that as in the case of alphanumerical database fields, indexes improve performances only when reading data using the SELECT command. In this case, you are making a number of updates on the table (INSERT, UPDATE, and DELETE); depending on the scenario, it could be less time consuming to drop and recreate the index after the updates.

You will import the Global_24h.csv file to the PostGIS database from NASA's Earth Observing System Data and Information System (EOSDIS).

You can copy the file from the dataset directory of the book for this chapter.

This file represents the active hotspots in the world detected by the Moderate Resolution Imaging Spectroradiometer (MODIS) satellites in the last 24 hours. For each row, there are the coordinates of the hotspot (latitude, longitude) in decimal degrees (in the WGS 84 spatial reference system, SRID = 4326), and a series of useful fields such as the acquisition date, acquisition time, and satellite type, just to name a few.

You will import only the active fire data scanned by the satellite type marked as T (Terra MODIS), and you will project it using the Spherical Mercator projection coordinate system (EPSG:3857; it is sometimes marked as EPSG:900913, where the number 900913 represents Google in 1337 speak, as it was first widely used by Google Maps).

The steps you need to follow to complete this recipe are as follows:

      $ cd ~/postgis_cookbook/data/chp01/
      $ head -n 5 Global_24h.csv

The output of the preceding command is as follows:

       <OGRVRTDataSource> 
         <OGRVRTLayer name="Global_24h"> 
         <SrcDataSource>Global_24h.csv</SrcDataSource> 
         <GeometryType>wkbPoint</GeometryType> 
         <LayerSRS>EPSG:4326</LayerSRS> 
           <GeometryField encoding="PointFromColumns"
            x="longitude" y="latitude"/> 
         </OGRVRTLayer> 
       </OGRVRTDataSource>

      $ ogrinfo global_24h.vrt Global_24h -fid 1

The output of the preceding command is as follows:

You can also try to open the virtual layer with a desktop GIS supporting a GDAL/OGR virtual driver such as Quantum GIS (QGIS). In the following screenshot, the Global_24h layer is displayed together with the shapefile of the countries that you can find in the dataset directory of the book:

The global_24h dataset over the countries layers and information of the selected features

      $ ogr2ogr -f PostgreSQL -t_srs EPSG:3857
      PG:"dbname='postgis_cookbook' user='me' password='mypassword'"
      -lco SCHEMA=chp01 global_24h.vrt -where "satellite='T'" 
      -lco GEOMETRY_NAME=the_geom

      $ pg_dump -t chp01.global_24h --schema-only -U me postgis_cookbook 
 
      CREATE TABLE global_24h ( 
        ogc_fid integer NOT NULL, 
        latitude character varying, 
        longitude character varying, 
        brightness character varying, 
        scan character varying, 
        track character varying, 
        acq_date character varying, 
        acq_time character varying, 
        satellite character varying, 
        confidence character varying, 
        version character varying, 
        bright_t31 character varying, 
        frp character varying, 
        the_geom public.geometry(Point,3857) 
      );

      postgis_cookbook=# SELECT f_geometry_column, coord_dimension,
      srid, type FROM geometry_columns 
      WHERE f_table_name = 'global_24h';

The output of the preceding command is as follows:

      postgis_cookbook=# SELECT count(*) FROM chp01.global_24h;

The output of the preceding command is as follows:

      postgis_cookbook=# SELECT ST_AsEWKT(the_geom)
      FROM chp01.global_24h LIMIT 1;

The output of the preceding command is as follows:

As mentioned in the GDAL documentation:

GDAL supports the reading and writing of nonspatial tabular data stored as a CSV file, but we need to use a virtual format to derive the geometry of the layers from attribute columns in the CSV file (the longitude and latitude coordinates for each point). For this purpose, you need to at least specify in the driver the path to the CSV file (the SrcDataSource element), the geometry type (the GeometryType element), the spatial reference definition for the layer (the LayerSRS element), and the way the driver can derive the geometric information (the GeometryField element).

There are many other options and reasons for using OGR virtual formats; if you are interested in developing a better understanding, please refer to the GDAL documentation available at http://www.gdal.org/drv_vrt.html.

After a virtual format is correctly created, the original flat nonspatial dataset is spatially supported by GDAL and software-based on GDAL. This is the reason why we can manipulate these files with GDAL commands such as ogrinfo and ogr2ogr, and with desktop GIS software such as QGIS.

Once we have verified that GDAL can correctly read the features from the virtual driver, we can easily import them in PostGIS using the popular ogr2ogr command-line utility. The ogr2ogr command has a plethora of options, so refer to its documentation at http://www.gdal.org/ogr2ogr.html for a more in-depth discussion.

In this recipe, you have just seen some of these options, such as:

The steps you need to follow to complete this recipe are as follows:

      $ ogr2ogr global_24h.shp global_24h.vrt

      $ shp2pgsql -G -I global_24h.shp
        chp01.global_24h_geographic > global_24h.sql

      $ head -n 20 global_24h.sql

The output of the preceding command is as follows:

      $ psql -U me -d postgis_cookbook -f global_24h.sql

      postgis_cookbook=# SELECT f_geography_column,   coord_dimension,
      srid, type FROM geography_columns   
      WHERE f_table_name = 'global_24h_geographic';

The output of the preceding command is as follows:

      $ ogrinfo PG:"dbname='postgis_cookbook' user='me'
        password='mypassword'" chp01.global_24h_geographic -fid 1

The output of the preceding command is as follows:

Now, open QGIS and try to add the new layer to the map. Navigate to Layer | Add Layer | Add PostGIS layers and provide the connection information, and then add the layer to the map as shown in the following screenshot:

The PostGIS command, shp2pgsql, allows the user to import a shapefile in the PostGIS database. Basically, it generates a PostgreSQL dump file that can be used to load data by running it from within PostgreSQL.

The SQL file will be generally composed of the following sections:

To get a complete list of the shp2pgsql command options and their meanings, just type the command name in the shell (or in the command prompt, if you are on Windows) and check the output.

There are GUI tools to manage data in and out of PostGIS, generally integrated into GIS desktop software such as QGIS. In the last chapter of this book, we will take a look at the most popular one.

The steps you need to follow to complete this recipe are as follows:

      $ ogr2ogr -f PostgreSQL -sql "SELECT ISO2, 
      NAME AS country_name FROM wborders WHERE REGION=2" -nlt 
      MULTIPOLYGON PG:"dbname='postgis_cookbook' user='me'
      password='mypassword'" -nln africa_countries 
      -lco SCHEMA=chp01 -lco GEOMETRY_NAME=the_geom wborders.shp

      postgis_cookbook=# SELECTST_AsText(the_geom) AS the_geom, bright_t31
FROM chp01.global_24h
ORDER BY bright_t31 DESC LIMIT 100;

The output of the preceding command is as follows:

      postgis_cookbook=# SELECT ST_AsText(f.the_geom) 
      AS the_geom, f.bright_t31, ac.iso2, ac.country_name
      FROM chp01.global_24h as f
      JOIN chp01.africa_countries as ac
      ON ST_Contains(ac.the_geom, ST_Transform(f.the_geom, 4326))
      ORDER BY f.bright_t31 DESCLIMIT 100;

The output of the preceding command is as follows:

You will now export the result of this query to a vector format supported by GDAL, such as GeoJSON, in the WGS 84 spatial reference using ogr2ogr:

      $ ogr2ogr -f GeoJSON -t_srs EPSG:4326 warmest_hs.geojson
      PG:"dbname='postgis_cookbook' user='me' password='mypassword'" -sql "
      SELECT f.the_geom as the_geom, f.bright_t31, 
             ac.iso2, ac.country_name
      FROM chp01.global_24h as f JOIN chp01.africa_countries as ac
      ON ST_Contains(ac.the_geom, ST_Transform(f.the_geom, 4326))
      ORDER BY f.bright_t31 DESC LIMIT 100"

      $ ogr2ogr -t_srs EPSG:4326 -f CSV -lco GEOMETRY=AS_XY 
      -lco SEPARATOR=TAB warmest_hs.csv PG:"dbname='postgis_cookbook' 
       user='me' password='mypassword'" -sql "
       SELECT f.the_geom, f.bright_t31,
         ac.iso2, ac.country_name 
       FROM chp01.global_24h as f JOIN chp01.africa_countries as ac 
       ON ST_Contains(ac.the_geom, ST_Transform(f.the_geom, 4326)) 
       ORDER BY f.bright_t31 DESC  LIMIT 100"

GDAL is an open source library that comes together with several command-line utilities, which let the user translate and process raster and vector geodatasets into a plethora of formats. In the case of vector datasets, there is a GDAL sublibrary for managing vector datasets named OGR (therefore, when talking about vector datasets in the context of GDAL, we can also use the expression OGR dataset).

When you are working with an OGR dataset, two of the most popular OGR commands are ogrinfo, which lists many kinds of information from an OGR dataset, and ogr2ogr, which converts the OGR dataset from one format to another.

It is possible to retrieve a list of the supported OGR vector formats using the -formats option on any OGR commands, for example, with ogr2ogr:

$ ogr2ogr --formats

The output of the preceding command is as follows:

Note that some formats are read-only, while others are read/write.

PostGIS is one of the supported read/write OGR formats, so it is possible to use the OGR API or any OGR commands (such as ogrinfo and ogr2ogr) to manipulate its datasets.

The ogr2ogr command has many options and parameters; in this recipe, you have seen some of the most notable ones such as -f to define the output format, -t_srs to reproject/transform the dataset, and -sql to define an (eventually spatial) query in the input OGR dataset.

When using ogrinfo and ogr2ogr together with the desired option and parameters, you have to define the datasets. When specifying a PostGIS dataset, you need a connection string that is defined as follows:

PG:"dbname='postgis_cookbook' user='me' password='mypassword'"

You can find more information about the ogrinfo and ogr2ogr commands on the GDAL website available at http://www.gdal.org.

If you need more information about the PostGIS driver, you should check its related documentation page available at http://www.gdal.org/drv_pg.html.

If you didn't follow all the other recipes, be sure to import the hotspots (Global_24h.csv) and the countries dataset (countries.shp) in PostGIS. The following is how to do it with ogr2ogr (you should import both the datasets in their original SRID, 4326, to make spatial operations faster):

      $ ogr2ogr -f PostgreSQL PG:"dbname='postgis_cookbook' 
      user='me' password='mypassword'" -lco SCHEMA=chp01 global_24h.vrt 
      -lco OVERWRITE=YES -lco GEOMETRY_NAME=the_geom -nln hotspots

      $ ogr2ogr -f PostgreSQL -sql "SELECT ISO2, NAME AS country_name 
      FROM wborders" -nlt MULTIPOLYGON PG:"dbname='postgis_cookbook' 
      user='me' password='mypassword'" -nln countries 
      -lco SCHEMA=chp01 -lco OVERWRITE=YES 
      -lco GEOMETRY_NAME=the_geom wborders.shp

The steps you need to follow to complete this recipe are as follows:

      postgis_cookbook=> SELECT c.country_name, MIN(c.iso2) 
      as iso2, count(*) as hs_count FROM chp01.hotspots as hs 
      JOIN chp01.countries as c ON ST_Contains(c.the_geom, hs.the_geom) 
      GROUP BY c.country_name ORDER BY c.country_name;

The output of the preceding command is as follows:

      $ ogr2ogr -f CSV hs_countries.csv 
      PG:"dbname='postgis_cookbook' user='me' password='mypassword'"
      -lco SCHEMA=chp01 -sql "SELECT c.country_name, MIN(c.iso2) as iso2, 
      count(*) as hs_count FROM chp01.hotspots as hs 
      JOIN chp01.countries as c ON ST_Contains(c.the_geom, hs.the_geom) 
      GROUP BY c.country_name ORDER BY c.country_name"
      postgis_cookbook=> COPY (SELECT c.country_name, MIN(c.iso2) as iso2, 
      count(*) as hs_count  FROM chp01.hotspots as hs
      JOIN chp01.countries as c  ON ST_Contains(c.the_geom, hs.the_geom)
      GROUP BY c.country_name  ORDER BY c.country_name) 
      TO '/tmp/hs_countries.csv' WITH CSV HEADER;

        #!/bin/bash 
        while IFS="," read country iso2 hs_count 
        do 
          echo "Generating shapefile $iso2.shp for country 
          $country ($iso2) containing $hs_count features." 
          ogr2ogr out_shapefiles/$iso2.shp
          PG:"dbname='postgis_cookbook' user='me' password='mypassword'"
          -lco SCHEMA=chp01 -sql "SELECT ST_Transform(hs.the_geom, 4326), 
          hs.acq_date, hs.acq_time, hs.bright_t31 
          FROM chp01.hotspots as hs JOIN chp01.countries as c 
          ON ST_Contains(c.the_geom, ST_Transform(hs.the_geom, 4326))  
          WHERE c.iso2 = '$iso2'" done < hs_countries.csv 

      chmod 775 export_shapefiles.sh
      mkdir out_shapefiles
      $ ./export_shapefiles.sh
      Generating shapefile AL.shp for country 
        Albania (AL) containing 66 features.
      Generating shapefile DZ.shp for country 
        Algeria (DZ) containing 361 features.
...
      Generating shapefile ZM.shp for country 
        Zambia (ZM) containing 1575 features.
      Generating shapefile ZW.shp for country 
        Zimbabwe (ZW) containing 179 features.

        @echo off 
        for /f "tokens=1-3 delims=, skip=1" %%a in (hs_countries.csv) do ( 
          echo "Generating shapefile %%b.shp for country %%a 
                (%%b) containing %%c features" 
          ogr2ogr .\out_shapefiles\%%b.shp 
          PG:"dbname='postgis_cookbook' user='me' password='mypassword'" 
          -lco SCHEMA=chp01 -sql "SELECT ST_Transform(hs.the_geom, 4326), 
          hs.acq_date, hs.acq_time, hs.bright_t31 
          FROM chp01.hotspots as hs JOIN chp01.countries as c 
          ON ST_Contains(c.the_geom, ST_Transform(hs.the_geom, 4326)) 
          WHERE c.iso2 = '%%b'" 
        ) 

      >mkdir out_shapefiles
      >export_shapefiles.bat
      "Generating shapefile AL.shp for country 
       Albania (AL) containing 66 features"
      "Generating shapefile DZ.shp for country 
       Algeria (DZ) containing 361 features"
      ...
      "Generating shapefile ZW.shp for country 
       Zimbabwe (ZW) containing 179 features"

      postgis_cookbook=# CREATE TABLE chp01.hs_uploaded
      (
        ogc_fid serial NOT NULL,
        acq_date character varying(80),
        acq_time character varying(80),
        bright_t31 character varying(80),
        iso2 character varying,
        upload_datetime character varying,
        shapefile character varying,
        the_geom geometry(POINT, 4326),
        CONSTRAINT hs_uploaded_pk PRIMARY KEY (ogc_fid)
      );

      $ brew install findutils

        #!/bin/bash 
        for f in `find out_shapefiles -name \*.shp -printf "%f\n"` 
        do 
          echo "Importing shapefile $f to chp01.hs_uploaded PostGIS
            table..." #, ${f%.*}" 
          ogr2ogr -append -update  -f PostgreSQL
          PG:"dbname='postgis_cookbook' user='me'
          password='mypassword'" out_shapefiles/$f 
          -nln chp01.hs_uploaded -sql "SELECT acq_date, acq_time,
          bright_t31, '${f%.*}' AS iso2, '`date`' AS upload_datetime,  
         'out_shapefiles/$f' as shapefile FROM ${f%.*}" 
        done 

      $ chmod 775 import_shapefiles.sh
      $ ./import_shapefiles.sh
      Importing shapefile DO.shp to chp01.hs_uploaded PostGIS table
      ...
      Importing shapefile ID.shp to chp01.hs_uploaded PostGIS table
      ...
      Importing shapefile AR.shp to chp01.hs_uploaded PostGIS table
      ......

Now, go to step 14.

        @echo off 
        for %%I in (out_shapefiles\*.shp*) do ( 
          echo Importing shapefile %%~nxI to chp01.hs_uploaded
          PostGIS table... 
 
          ogr2ogr -append -update  -f PostgreSQL
          PG:"dbname='postgis_cookbook' user='me'
          password='password'" out_shapefiles/%%~nxI 

          -nln chp01.hs_uploaded -sql "SELECT acq_date, acq_time, 
          bright_t31, '%%~nI' AS iso2, '%date%' AS upload_datetime, 
          'out_shapefiles/%%~nxI' as shapefile FROM %%~nI" 
        ) 

      >import_shapefiles.bat
      Importing shapefile AL.shp to chp01.hs_uploaded PostGIS table...
      Importing shapefile AO.shp to chp01.hs_uploaded PostGIS table...
      Importing shapefile AR.shp to chp01.hs_uploaded PostGIS table......

      postgis_cookbook=# SELECT upload_datetime,
      shapefile, ST_AsText(wkb_geometry)
      FROM chp01.hs_uploaded WHERE ISO2='AT';

The output of the preceding command is as follows:

      $ ogrinfo PG:"dbname='postgis_cookbook' user='me'
      password='mypassword'"
      chp01.hs_uploaded -where "iso2='AT'"

The output of the preceding command is as follows:

You could implement both the data flows (processing shapefiles out from PostGIS, and then into it again) thanks to the power of the ogr2ogr GDAL command.

You have been using this command in different forms and with the most important input parameters in other recipes, so you should now have a good understanding of it.

Here, it is worth mentioning the way OGR lets you export the information related to the current datetime and the original shapefile name to the PostGIS table. Inside the import_shapefiles.sh (Linux, OS X) or the import_shapefiles.bat (Windows) scripts, the core is the line with the ogr2ogr command (here is the Linux version):

ogr2ogr -append -update  -f PostgreSQL PG:"dbname='postgis_cookbook' user='me' password='mypassword'" out_shapefiles/$f -nln chp01.hs_uploaded -sql "SELECT acq_date, acq_time, bright_t31, '${f%.*}' AS iso2, '`date`' AS upload_datetime, 'out_shapefiles/$f' as shapefile FROM ${f%.*}"

Thanks to the -sql option, you can specify the two additional fields, getting their values from the system date command and the filename that is being iterated from the script.

The steps you need to follow to complete this recipe are as follows:

      $ shp2pgsql -I -d -s 4326 -W LATIN1 -g the_geom countries.shp
      chp01.countries > countries.sql
      $ psql -U me -d postgis_cookbook -f countries.sql

      $ ogr2ogr -f PostgreSQL PG:"dbname='postgis_cookbook' user='me'
      password='mypassword'"
      -lco SCHEMA=chp01 countries.shp -nlt MULTIPOLYGON -lco OVERWRITE=YES
      -lco GEOMETRY_NAME=the_geom

      postgis_cookbook=> SELECT subregion,
        ST_Union(the_geom) AS the_geom, SUM(pop2005) AS pop2005
        FROM chp01.countries GROUP BY subregion;

      $ pgsql2shp -f subregions.shp -h localhost -u me -P mypassword
      postgis_cookbook "SELECT MIN(subregion) AS subregion, 
      ST_Union(the_geom) AS the_geom, SUM(pop2005) AS pop2005 
      FROM chp01.countries GROUP BY subregion;" 
 
      Initializing... 
      Done (postgis major version: 2). 
      Output shape: Polygon 
      Dumping: X [23 rows].

Visualization in QGIS of the classification of subregions based on population and information of the selected feature

You have exported the results of a spatial query to a shapefile using the pgsql2shp PostGIS command. The spatial query you have used aggregates fields using the SUM PostgreSQL function for summing country populations in the same subregion, and the ST_Union PostGIS function to aggregate the corresponding geometries as a geometric union.

The pgsql2shp command allows you to export PostGIS tables and queries to shapefiles. The options you need to specify are quite similar to the ones you use to connect to PostgreSQL with psql. To get a full list of these options, just type pgsql2shp in your command prompt and read the output.

We need the following in place before we can proceed with the steps required for the recipe:

      $ sudo apt-get install osm2pgsql

      $ osm2pgsqlosm2pgsql SVN version 0.80.0 (32bit id space)

      postgres=# CREATE DATABASE rome OWNER me;
      postgres=# \connect rome;
      rome=# create extension postgis;

      $ sudo apt-get update
      $ sudo apt-get install postgresql-contrib-9.6
      $ psql -U me -d romerome=# CREATE EXTENSION hstore;

The steps you need to follow to complete this recipe are as follows:

      $ osm2pgsql -d rome -U me --hstore map.osm
      osm2pgsql SVN version 0.80.0 (32bit id space)Using projection
      SRS 900913 (Spherical Mercator)Setting up table: 
      planet_osm_point...All indexes on planet_osm_polygon created 
      in 1sCompleted planet_osm_polygonOsm2pgsql took 3s overall

      rome=# SELECT f_table_name, f_geometry_column,
      coord_dimension, srid, type FROM geometry_columns;

The output of the preceding command is shown here:

      rome=# CREATE VIEW rome_trees AS SELECT way, tags 
      FROM planet_osm_polygon WHERE (tags -> 'landcover') = 'trees';

OpenStreetMap is a popular collaborative project for creating a free map of the world. Every user participating in the project can edit data; at the same time, it is possible for everyone to download those datasets in .osm datafiles (an XML format) under the terms of the Open Data Commons Open Database License (ODbL) at the time of writing.

The osm2pgsql command is a command-line tool that can import .osm datafiles (eventually zipped) to the PostGIS database. To use the command, it is enough to give the PostgreSQL connection parameters and the .osm file to import.

It is possible to import only features that have certain tags in the spatial database, as defined in the default.style configuration file. You can decide to comment in or out the OSM tagged features that you would like to import, or not, from this file. The command by default exports all the nodes and ways to linestring, point, and geometry PostGIS geometries.

It is highly recommended to enable hstore support in the PostgreSQL database and use the -hstore option of osm2pgsql when importing the data. Having enabled this support, the OSM tags for each feature will be stored in a hstore PostgreSQL data type, which is optimized for storing (and retrieving) sets of key/values pairs in a single field. This way, it will be possible to query the database as follows:

SELECT way, tags FROM planet_osm_polygon WHERE (tags -> 'landcover') = 'trees';

We need the following in place before we can proceed with the steps required for the recipe:

      $ shp2pgsql -I -d -s 4326 -W LATIN1 -g the_geom countries.shp
      chp01.countries > countries.sql
      $ psql -U me -d postgis_cookbook -f countries.sql

The steps you need to follow to complete this recipe are as follows:

      $ gdalinfo worldclim/tmax09.bil

      Driver: EHdr/ESRI .hdr Labelled
      Files: worldclim/tmax9.bil
worldclim/tmax9.hdr
      Size is 2160, 900
      Coordinate System is:
      GEOGCS[""WGS 84"",
        DATUM[""WGS_1984"",
          SPHEROID[""WGS 84"",6378137,298.257223563,
            AUTHORITY[""EPSG"",""7030""]],
          TOWGS84[0,0,0,0,0,0,0],
          AUTHORITY[""EPSG"",""6326""]],
        PRIMEM[""Greenwich"",0,
          AUTHORITY[""EPSG"",""8901""]],
        UNIT[""degree"",0.0174532925199433,
        AUTHORITY[""EPSG"",""9108""]],
      AUTHORITY[""EPSG"",""4326""]]
      Origin = (-180.000000000000057,90.000000000000000)
      Pixel Size = (0.166666666666667,-0.166666666666667)
      Corner Coordinates:
Upper Left (-180.0000000, 90.0000000) (180d 0'' 0.00""W, 90d
                                               0'' 0.00""N)
Lower Left (-180.0000000, -60.0000000) (180d 0'' 0.00""W, 60d
                                                0'' 0.00""S)
Upper Right ( 180.0000000, 90.0000000) (180d 0'' 0.00""E, 90d
                                                0'' 0.00""N)
Lower Right ( 180.0000000, -60.0000000) (180d 0'' 0.00""E, 60d 
                                                 0'' 0.00""S)
Center ( 0.0000000, 15.0000000) ( 0d 0'' 0.00""E, 15d
                                          0'' 0.00""N)
        Band 1 Block=2160x1 Type=Int16, ColorInterp=Undefined
        Min=-153.000 Max=441.000
        NoData Value=-9999

      $ raster2pgsql -I -C -F -t 100x100 -s 4326
      worldclim/tmax01.bil chp01.tmax01 > tmax01.sql
      $ psql -d postgis_cookbook -U me -f tmax01.sql

If you are in Linux, you may pipe the two commands in a unique line:

      $ raster2pgsql -I -C -M -F -t 100x100 worldclim/tmax01.bil 
      chp01.tmax01 | psql -d postgis_cookbook -U me -f tmax01.sql

      $ pg_dump -t chp01.tmax01 --schema-only -U me postgis_cookbook
      ...
      CREATE TABLE tmax01 (
        rid integer NOT NULL,
        rast public.raster,
        filename text,
        CONSTRAINT enforce_height_rast CHECK (
          (public.st_height(rast) = 100)
        ),
        CONSTRAINT enforce_max_extent_rast CHECK (public.st_coveredby
          (public.st_convexhull(rast), ''0103...''::public.geometry)
        ),
        CONSTRAINT enforce_nodata_values_rast CHECK (
          ((public._raster_constraint_nodata_values(rast)
            )::numeric(16,10)[] = ''{0}''::numeric(16,10)[])
          ),
        CONSTRAINT enforce_num_bands_rast CHECK (
          (public.st_numbands(rast) = 1)
        ),
        CONSTRAINT enforce_out_db_rast CHECK (
          (public._raster_constraint_out_db(rast) = ''{f}''::boolean[])
          ),
        CONSTRAINT enforce_pixel_types_rast CHECK (
          (public._raster_constraint_pixel_types(rast) = 
           ''{16BUI}''::text[])
          ),
        CONSTRAINT enforce_same_alignment_rast CHECK (
          (public.st_samealignment(rast, ''01000...''::public.raster)
        ),
        CONSTRAINT enforce_scalex_rast CHECK (
          ((public.st_scalex(rast))::numeric(16,10) = 
            0.166666666666667::numeric(16,10))
           ),
        CONSTRAINT enforce_scaley_rast CHECK (
          ((public.st_scaley(rast))::numeric(16,10) = 
            (-0.166666666666667)::numeric(16,10))
          ),
        CONSTRAINT enforce_srid_rast CHECK ((public.st_srid(rast) = 0)),
        CONSTRAINT enforce_width_rast CHECK ((public.st_width(rast) = 100))
      );

      postgis_cookbook=# SELECT * FROM raster_columns;

      postgis_cookbook=# SELECT count(*) FROM chp01.tmax01;

The output of the preceding command is as follows:

      count
      -------
      198
      (1 row)

      gdalinfo PG":host=localhost port=5432 dbname=postgis_cookbook
      user=me password=mypassword schema='chp01' table='tmax01'"

      gdalinfo PG":host=localhost port=5432 dbname=postgis_cookbook
      user=me password=mypassword schema='chp01' table='tmax01' mode=2"

      $ ogr2ogr temp_grid.shp PG:"host=localhost port=5432 
      dbname='postgis_cookbook' user='me' password='mypassword'" 
      -sql "SELECT rid, filename, ST_Envelope(rast) as the_geom 
      FROM chp01.tmax01"

      SELECT * FROM (
        SELECT c.name, ST_Value(t.rast,
          ST_Centroid(c.the_geom))/10 as tmax_jan FROM chp01.tmax01 AS t 
        JOIN chp01.countries AS c 
        ON ST_Intersects(t.rast, ST_Centroid(c.the_geom))
      ) AS foo 
      ORDER BY tmax_jan LIMIT 10;

The output is as follows:

The raster2pgsql command is able to load any raster formats supported by GDAL in PostGIS. You can have a format list supported by your GDAL installation by typing the following command:

$ gdalinfo --formats

In this recipe, you have been importing one raster file using some of the most common raster2pgsql options:

$ raster2pgsql -I -C -F -t 100x100 -s 4326 worldclim/tmax01.bil chp01.tmax01 > tmax01.sql

The -I option creates a GIST spatial index for the raster column. The -C option will create the standard set of constraints after the rasters have been loaded. The -F option will add a column with the filename of the raster that has been loaded. This is useful when you are appending many raster files to the same PostGIS raster table. The -s option sets the raster's SRID.

If you decide to include the -t option, then you will cut the original raster into tiles, each inserted as a single row in the raster table. In this case, you decided to cut the raster into 100 x 100 tiles, resulting in 198 table rows in the raster table.

Another important option is -R, which will register the raster as out-of-db; in such a case, only the metadata will be inserted in the database, while the raster will be out of the database.

The raster table contains an identifier for each row, the raster itself (eventually one of its tiles, if using the -t option), and eventually the original filename, if you used the -F option, as in this case.

You can analyze the PostGIS raster using SQL commands or the gdalinfo command. Using SQL, you can query the raster_columns view to get the most significant raster metadata (spatial reference, band number, scale, block size, and so on).

With gdalinfo, you can access the same information, using a connection string with the following syntax:

gdalinfo PG":host=localhost port=5432 dbname=postgis_cookbook user=me password=mypassword schema='chp01' table='tmax01' mode=2"

The mode parameter is not influential if you loaded the whole raster as a single block (for example, if you did not specify the -t option). But, as in the use case of this recipe, if you split it into tiles, gdalinfo will see each tile as a single subdataset with the default behavior (mode=1). If you want GDAL to consider the raster table as a unique raster dataset, you have to specify the mode option and explicitly set it to 2.

Be sure to have all the original raster datasets you have been using for the previous recipe.

The steps you need to follow to complete this recipe are as follows:

      $ raster2pgsql -d -I -C -M -F -t 100x100 -s 4326 
      worldclim/tmax*.bil chp01.tmax_2012 > tmax_2012.sql
      $ psql -d postgis_cookbook -U me -f tmax_2012.sql

      postgis_cookbook=# SELECT r_raster_column, srid,
      ROUND(scale_x::numeric, 2) AS scale_x, 
      ROUND(scale_y::numeric, 2) AS scale_y, blocksize_x, 
      blocksize_y, num_bands, pixel_types, nodata_values, out_db 
      FROM raster_columns where r_table_schema='chp01' 
      AND r_table_name ='tmax_2012';

      SELECT rid, (foo.md).* 
      FROM (SELECT rid, ST_MetaData(rast) As md 
      FROM chp01.tmax_2012) As foo;

The output of the preceding command is as shown here:

      postgis_cookbook=# SELECT COUNT(*) AS num_raster, 
      MIN(filename) as original_file FROM chp01.tmax_2012
GROUP BY filename ORDER BY filename;

      SELECT REPLACE(REPLACE(filename, 'tmax', ''), '.bil', '') AS month, 
      (ST_VALUE(rast, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) AS tmax 
      FROM chp01.tmax_2012 
      WHERE rid IN (
        SELECT rid FROM chp01.tmax_2012 
        WHERE ST_Intersects(ST_Envelope(rast),
              ST_SetSRID(ST_Point(12.49, 41.88), 4326))
      )
      ORDER BY month;

The output of the preceding command is as shown here:

      $ gdalbuildvrt -separate tmax_2012.vrt worldclim/tmax*.bil

      <VRTDataset rasterXSize="2160" rasterYSize="900">
        <SRS>GEOGCS...</SRS>
        <GeoTransform> 
          -1.8000000000000006e+02, 1.6666666666666699e-01, ...
        </GeoTransform>
        <VRTRasterBand dataType="Int16" band="1">
          <NoDataValue>-9.99900000000000E+03</NoDataValue>
          <ComplexSource>
            <SourceFilename relativeToVRT="1">
              worldclim/tmax01.bil
            </SourceFilename>
            <SourceBand>1</SourceBand>
            <SourceProperties RasterXSize="2160" RasterYSize="900"
             DataType="Int16" BlockXSize="2160" BlockYSize="1" />
            <SrcRect xOff="0" yOff="0" xSize="2160" ySize="900" />
            <DstRect xOff="0" yOff="0" xSize="2160" ySize="900" />
            <NODATA>-9999</NODATA>
          </ComplexSource>
        </VRTRasterBand>
      <VRTRasterBand dataType="Int16" band="2">
      ...

      $ gdalinfo tmax_2012.vrt

The output of the preceding command is as follows:

      $ raster2pgsql -d -I -C -M -F -t 100x100 -s 4326 tmax_2012.vrt 
      chp01.tmax_2012_multi > tmax_2012_multi.sql
      $ psql -d postgis_cookbook -U me -f tmax_2012_multi.sql

     postgis_cookbook=# SELECT r_raster_column, srid, blocksize_x,
     blocksize_y, num_bands, pixel_types 
     from raster_columns where r_table_schema='chp01' 
     AND r_table_name ='tmax_2012_multi';

      postgis_cookbook=# SELECT
      (ST_VALUE(rast, 1, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS jan,
      (ST_VALUE(rast, 2, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS feb,
      (ST_VALUE(rast, 3, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10)
      AS mar,
      (ST_VALUE(rast, 4, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS apr,
      (ST_VALUE(rast, 5, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS may,
      (ST_VALUE(rast, 6, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS jun,
      (ST_VALUE(rast, 7, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS jul,
      (ST_VALUE(rast, 8, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS aug,
      (ST_VALUE(rast, 9, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS sep,
      (ST_VALUE(rast, 10, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS oct,
      (ST_VALUE(rast, 11, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS nov,
      (ST_VALUE(rast, 12, ST_SetSRID(ST_Point(12.49, 41.88), 4326))/10) 
      AS dec 
      FROM chp01.tmax_2012_multi WHERE rid IN (
        SELECT rid FROM chp01.tmax_2012_multi
        WHERE ST_Intersects(rast, ST_SetSRID(ST_Point(12.49, 41.88), 4326))
      );

The output of the preceding command is as follows:

You can import raster datasets in PostGIS using the raster2pgsql command.

In a scenario where you have multiple rasters representing the same variable at different times, as in this recipe, it makes sense to store all of the original rasters as a single table in PostGIS. In this recipe, we have the same variable (average maximum temperature) represented by a single raster for each month. You have seen that you could proceed in two different ways:

You need the following in place before you can proceed with the steps required for the recipe:

      $ gdalinfo --formats | grep -i postgis

The output of the preceding command is as follows:

      PostGISRaster (rw): PostGIS Raster driver

      $ gdalinfo PG:"host=localhost port=5432
      dbname='postgis_cookbook' user='me' password='mypassword'
      schema='chp01' table='tmax_2012_multi' mode='2'"

The steps you need to follow to complete this recipe are as follows:

      $ gdal_translate -b 1 -b 2 -b 3 -b 4 -b 5 -b 6
      PG:"host=localhost port=5432 dbname='postgis_cookbook'
      user='me' password='mypassword' schema='chp01'
      table='tmax_2012_multi' mode='2'" tmax_2012_multi_123456.tif

      postgis_cookbook=# SELECT ST_Extent(the_geom) 
      FROM chp01.countries WHERE name = 'Italy';

The output of the preceding command is as follows:

      $ gdal_translate -projwin 6.619 47.095 18.515 36.649
      PG:"host=localhost port=5432 dbname='postgis_cookbook'
      user='me' password='mypassword' schema='chp01'
      table='tmax_2012_multi' mode='2'" tmax_2012_multi.tif

      gdalwarp -t_srs EPSG:3857 PG:"host=localhost port=5432 
      dbname='postgis_cookbook' user='me' password='mypassword' 
      schema='chp01' table='tmax_2012_multi' mode='2'" 
      tmax_2012_multi_3857.tif

Both gdal_translate and gdalwarp can transform rasters from a PostGIS raster to all GDAL-supported formats. To get a complete list of the supported formats, you can use the --formats option of GDAL's command line as follows:

$ gdalinfo --formats

For both these GDAL commands, the default output format is GeoTiff; if you need a different format, you must use the -of option and assign to it one of the outputs produced by the previous command line.

In this recipe, you have tried some of the most common options for these two commands. As they are complex tools, you may try some more command options as a bonus step.

To get a better understanding, you should check out the excellent documentation on the GDAL website: