Outputting rasters using GDAL

We have seen ogr2ogr in action already, and now it is time to give two GDAL tools a spin: gdal_translate and gdalwarp. The former is a translation utility that can change the format of a raster datasource while the latter is a reprojection utility. Therefore, in order to extract a reprojected raster from a PostGIS database, we need to perform two steps:

  1. Extract the raster using gdal_translate
  2. Perform a reprojection using gdalwarp

In order to get some help with both utilities, simply type the following in the console:

gdal_translate

You could also use the following:

gdalwarp

The most important parameters of gdal_translate for our scenario are:

  • -of: Specifies the output format (use gdal_translate -formats to obtain information on the supported formats).
  • -outsize xsize[%] ysize[%]: Specifies the output raster size. It can be expressed in pixels or a percentage. If it is not used 100%, then the original size is assumed
  • PG (PgSQL connection opts): PgSQL connection options. It is possible to pass a more complex where clause using a command line where you switch to ogrinfo.
  • where: Filtering using a where clause.
  • -projwin upper left X, upper left Y, lower right X, lower right Y: Provides a means for clipping the output raster using projected coordinates. If for some reason you prefer using pixel coordinates, use the -srcwin parameter.

So, let's export our Natural Earth raster to GeoTiff first, making sure that we export the top-right section and that we shrink it by half on the way:

gdal_translate PG:"host=localhost port=5434 user=postgres dbname=mastering_postgis schema=data_import table=gray_50m_partial where='filename='gray_50m_partial_bl.tif'' mode=2" -of GTiff -outsize 50% 50% gray_50m_partial_small.tiff

You should see similar output to the following:

Input file size is 5400, 2700
0...10...20...30...40...50...60...70...80...90...100 - done.
The following is taken from the PostGIS documentation: The mode=2 is required for tiled rasters and was added in PostGIS 2.0 and GDAL 1.8 drivers. This does not exist in GDAL 1.7 drivers.

The where clause can also contain spatial filtering. The following command exports a tile that intersects with some declared bounds - in this case, it's roughly Poland's bounding box:

gdal_translate PG:"host=localhost port=5434 user=postgres dbname=mastering_postgis schema=data_import table=gray_50m_partial where='ST_Intersects(rast, ST_MakeEnvelope(14,49,24,55,4326))' mode=2" -of GTiff -outsize 50% 50% gray_50m_partial_small.tiff

Although this example uses an ST_Intersects function, it does not clip the raster to the bounds used. Instead, it just picks a tile that intersects them. However, we still have some options to clip a raster, using, for example, the aforementioned -projwin parameter. To expand on the previous example, let's now crop the territory of Poland out of the larger raster:

gdal_translate PG:"host=localhost port=5434 user=postgres dbname=mastering_postgis schema=data_import table=gray_50m_partial where='ST_Intersects(rast, ST_MakeEnvelope(14,49,24,55,4326))' mode=2" -of GTiff -projwin 14 55 24 49 poland.tiff
Should you require some more sophisticated filtering or better control over your output, you can use a view or a temporary table to prepare the data. Unfortunately, with GDAL we cannot issue an SQL select query as we were able to do with ogr2ogr.

Having exported a raster to a file, changing its projection is now a formality. Let's reproject our Poland raster to EPSG:2180 - Polish Coordinate System 92:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:2180 poland.tiff poland_2180.tiff

Should you need to adjust the resampling method when using gdalwarp, you can do so using the -r parameter. You could also use gdalwarp to crop a raster; use a -te xmin ymin xmax ymax to achieve this.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.
Reset