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:
- Extract the raster using gdal_translate
- 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 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
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.