This recipe will guide you through the importing 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, with one band for each of the original rasters, and then load the multiband raster to a raster table. To accomplish this, you will use the GDAL gdalbuildvrt command and then load the data to PostGIS with raster2pgsql.
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:
- Import all the maximum average temperature rasters in a single PostGIS raster table using
raster2pgsqland thenpsql(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
- Check how the table was created in PostGIS, querying the
raster_columnstable. 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';
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:
If you now query the table, you would be able to derive the month for each raster row only from the
original_filecolumn. In the table, you have imported 198 distinct records (rasters) for each of the 12 original files (we divided them into 100 x 100 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;
- With this approach, using the
filenamefield, you could use theST_ValuePostGIS 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;The output of the preceding command is as shown here:
- A different approach is to store each month value in a different raster band. The
raster2pgsqlcommand doesn't let you load to different bands in an existing table. But, you can use GDAL by combining thegdalbuildvrtand thegdal_translatecommands. First, usegdalbuildvrtto create a new virtual raster composed of 12 bands, one for each month:
$ gdalbuildvrt -separate tmax_2012.vrt worldclim/tmax*.bil- Analyze the
tmax_2012.vrtXML 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">
...- Now, with
gdalinfo, analyze this output virtual raster to check if it is effectively composed of 12 bands:
$ gdalinfo tmax_2012.vrtThe output of the preceding command is as follows:
- 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
raster2pgsqlcommand:
$ 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
- Query the
raster_columnsview to get some indicators for the imported raster. Note that thenum_bandsis 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';
- Now, let's try to produce the same output as the query using 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:
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 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:
- 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 raster2pgsqloption). - 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 theraster2pgsqlcommand.