Importing multiple rasters at a time

This recipe will guide you through the import of multiple rasters at a time.

You will first import some different single band rasters to a unique single band raster table using the raster2pgsql command.

Then, you will try an alternative approach, merging the original single band rasters in a virtual raster, having one band for each of the original rasters, and then load the multiband raster to a raster table. For accomplishing this, you will use the GDAL gdalbuildvrt command and then load the data to PostGIS with raster2pgsql.

Getting ready

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

How to do it...

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

1. Import all the maximum average temperature rasters in a single PostGIS raster table using raster2pgsql and then psql (eventually, pipe the two commands if you are in Linux):

   $ 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
   

2. Check how the table was created in PostGIS, querying the raster_columns table. Here we are querying only some significant fields:

   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';
    r_raster_column | srid | scale_x | scale_y | blocksize_x | blocksize_y | num_bands | pixel_types | nodata_values | out_db-----------------+------+---------+---------+-------------+-------------+-----------+-------------+---------------+--------rast            | 4326 |    0.17 |   -0.17 |         100 |         100 |         1 | {16BUI}     | {0}           | {f}
   (1 row)
   

3. Check some raster statistics using the ST_MetaData function:

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

   Note

   Note that there is a different metadata for each raster record loaded in the table.

4. If you now query the table, you would be able to derive the month for each raster row only from the original_file column. In the table, you have imported 198 distinct records (raster) for each of the 12 original files (we divided them into 100x100 blocks, if you remember). Test this with the following query:

   postgis_cookbook=# SELECT COUNT(*) AS num_raster, MIN(filename) as original_file FROM chp01.tmax_2012
   GROUP BY filename ORDER BY filename;
    num_raster | original_file
   ------------+---------------
           198 | tmax01.bil
           198 | tmax02.bil
           198 | tmax03.bil
           198 | tmax04.bil
           198 | tmax05.bil
           198 | tmax06.bil
           198 | tmax07.bil
           198 | tmax08.bil
           198 | tmax09.bil
           198 | tmax10.bil
           198 | tmax11.bil
           198 | tmax12.bil
   (12 rows)
   

5. With this approach, using the filename field, you could use the ST_Value PostGIS raster function to get the average monthly maximum temperature of a certain geographic zone for the whole year:

   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;
    month | tmax
   -------+------
    01    | 11.8
    02    | 13.2
    03    | 15.3
    04    | 18.5
    05    | 22.9
    06    | 27
    07    | 30
    08    | 29.8
    09    | 26.4
    10    | 21.7
    11    | 16.6
    12    | 12.9
   (12 rows)
   

6. A different approach is to store each month value in a different raster band. The raster2pgsql command doesn't let you load to different bands in an existing table. But, you can use GDAL by combining the gdalbuildvrt and the gdal_translate commands. First, use gdalbuildvrt to create a new virtual raster composed of 12 bands, one for each month:

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

7. Analyze the tmax_2012.vrt xml file with a text editor. It should have a virtual band (VRTRasterBand) for each physical raster pointing to it:

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

8. Now, with gdalinfo, analyze this output virtual raster to check if it is effectively composed of 12 bands:

   $ gdalinfo tmax_2012.vrt
   

   The output of the preceding command is as follows:

   Driver: VRT/Virtual Raster
   Files: tmax_2012.vrt
          worldclim/tmax01.bil
          worldclim/tmax02.bil
          ...
          worldclim/tmax12.bil
   Size is 2160, 900
   Coordinate System is:
   GEOGCS["WGS 84",
       DATUM["WGS_1984",
           ...
   Band 1 Block=128x128 Type=Int16, ColorInterp=Undefined
     Min=-478.000 Max=418.000
     NoData Value=-9999
   Band 2 Block=128x128 Type=Int16, ColorInterp=Undefined
     Min=-421.000 Max=414.000
     NoData Value=-9999
   ...
   

9. Import the virtual raster composed of 12 bands, each referring to one of the 12 original rasters, to a PostGIS raster table composed of 12 bands. For this purpose, you can use the raster2pgsql command:

   $ 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
   

10. Query the raster_columns view to get some indicators for the imported raster. Note as the num_bands is now 12:

        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';
     r_raster_column | srid | blocksize_x | blocksize_y | num_bands |                                pixel_types
    -----------------+------+-------------+-------------+-----------+---------------------------------------------------------------------------
     rast            | 4326 |         100 |         100 |        12 | {16BSI,16BSI,16BSI,16BSI,16BSI,16BSI,16BSI,16BSI,16BSI,16BSI,16BSI,16BSI}
    

11. Now, let's try to produce the same output of the query with the previous approach. This time, given the table structure, we keep the results in a single row:

        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:

    jan | feb | mar | apr | may | jun | jul | aug | sep | oct | nov | dec
    ----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+----
    11.8| 13.2| 15.3| 18.5| 22.9| 27  | 30  | 29.8| 26.4| 21.7| 16.6|12.9
    (1 row)
    

How it works...

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

Note

The GDAL PostGIS raster so far does not support writing operations; therefore, for now, you cannot use GDAL commands such as gdal_translate and gdalwarp.

This is going to change in the near future, so you may have such an extra option when you are reading this chapter.

In a scenario where you have multiple rasters representing the same variable at a time, 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:

1. Append each single raster (representing a different month) to the same PostGIS single band raster table and derive the information related to the month from the value in the filename column (added to the table using the -F raster2pgsql option).

2. Generate a multiband raster using gdalbuildvrt (one raster with 12 bands, one for each month), and import it in a single multiband PostGIS table using the raster2pgsql command.