From the worldclim.org website, download the current raster data (http://www.worldclim.org/current) for min and max temperatures (only the raster for max temperatures will be used for this recipe). Alternatively, use the ones provided in the book datasets (data/chp01). Each of the two archives (data/tmax_10m_bil.zip and data/tmin_10m_bil.zip) contain 12 rasters in the BIL format, one for each month. You can look for more information at http://www.worldclim.org/formats.

Extract the two archives to a directory named worldclim in your working directory.

Rename each raster dataset to a name format having two digits for the month, for example, tmax1.bil and tmax1.hdr will become tmax01.bil and tmax01.hdr.

If you still haven't loaded the countries shapefile to PostGIS from a previous recipe, do it using the ogr2ogr or shp2pgsql commands. The following is the shp2pgsql syntax:

$ 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

Get information about one of the rasters using the gdalinfo command-line tool 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, 90d0' 0.00"N)
  Lower Left  (-180.0000000, -60.0000000) (180d 0' 0.00"W, 60d0' 0.00"S)
  Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d0' 0.00"N)
  Lower Right ( 180.0000000, -60.0000000) (180d 0' 0.00"E, 60d0' 0.00"S)
  Center      (   0.0000000,  15.0000000) (  0d 0' 0.00"E, 15d0' 0.00"N)
Band 1 Block=2160x1 Type=Int16, ColorInterp=Undefined
Min=-153.000 Max=441.000
NoData Value=-9999

The gdalinfo command provides a series of useful information about the raster, for example, the GDAL driver being used to read it, the files composing it (in this case, two files with a .bil and .hdr extension), the size in pixels (2160 x 900), the spatial reference (WGS 84), the geographic extents, the origin and the pixel size (needed to correctly georeference the raster), and for each raster band (just one in the case of this file), some statistical information like the min and max value (-153.000 and 441.000, corresponding to a temperature of -15.3 °C and 44.1 °C. Values are expressed as temperature * 10 in °C, according to the documentation available at worldclim.org).

Use the raster2pgsql file to generate the .sql dump file and then import the raster in PostGIS:

$ 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/tmax1.bil chp01.tmax01 | psql -d postgis_cookbook -U me -f tmax01.sql

Check how the new table has been created in PostGIS:

$ 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))
);

Check if a record for this PostGIS raster appears in the raster_columns metadata view, and note the main metadata information that has been stored there such as schema, name, raster column name (default is raster), SRID, scale (for x and y), block size (for x and y), band numbers (1), band types (16BUI), zero data values (0), and db storage type (out_db is false, as we have stored the raster bytes in the database; you could have used the -R option to register the raster as an out-of-db filesystem):

postgis_cookbook=# SELECT * FROM raster_columns;

If you have followed this recipe from the beginning, you should now have 198 rows in the raster table, with each row representing one raster block size (100 x 100 pixels blocks, as indicated with the -t raster2pgsql option):

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

The output of the preceding command is as follows:

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

Try to open the raster table with gdalinfo. You should see the same information you got from gdalinfo when you were analyzing the original BIL file. The only difference is the block size, as you moved to a smaller one (100x100) from the original (2160x900). That's why the original file has been split into several datasets (198):

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

The gdalinfo command reads the PostGIS raster as being composed of multiple raster subdatasets (198, one for each row in the table). You still have the possibility of reading the whole table as a single raster, using the mode=2 option in the PostGIS raster connection string (mode=1 is the default). Check the difference:

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

You can easily obtain a visual representation of those blocks by converting the extent of all the 198 rows in the tmax01 table (each representing a raster block) to a shapefile using ogr2ogr:

$ 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"

Now, try to open the raster table with QGIS (at the time of writing, one of the few Desktop GIS tools had support for it) together with the blocks shapefile generated in the previous steps (temp_grid.shp). You should see something like the following screenshot:

If you are using QGIS, you need the PostGIS raster QGIS plugin to read and write the PostGIS raster. This plugin makes it possible to add single rows of the tmax01 table as a single raster or the whole table. Try to add a couple of rows.

As the last bonus step, you will select the 10 countries with the lowest average max temperature in January (using the centroid of the polygon representing the country):

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:

name                                        | tmax_jan
--------------------------------------------+----------
Greenland                                   |    -29.8
 ...
Korea                                       |     -8.5
Democratic People's Republic of Kyrgyzstan  |     -7.9
Finland                                     |     -6.8
(10 rows)