We had this big raster that we needed to chop up into tiles and only extract a portion of for load into PostGIS. raster2pgsql doesn't currently have an option to pull just a portion of a raster and also we don't have the windows raster2pgsql compiled with MrSID support. Luckily
GDAL commandline gdal_translate has this switch that allows you to specify a chunk of a raster to grab based on a projected or unprojected box window.
We wanted to grab just that portion that is part of boston and chunk it into bite size pieces. What we needed was a grid generator similar to
what we described a while back in Map Dicing and other stuff
that would chop our neighborhood into bite sized tiles we could then use to generate the relevant gdal_translate command.
Instead of using temp tables and all that stuff, we decided to try with the ST_Tile raster function. Creating an empty raster and then tiling it.
Note the repurposing: Creating a raster with no bands to accomplish a task that has nothing to do with rasters, so that we can then apply it to something to do with rasters. Gridding is a surprisingly common step in a lot of spatial processing workflows.
Here is the SQL to do it and we'll explain in a separate article in more detail.
WITH dims As (
SELECT 10 As perpixel
, 1024 AS tile_width
, 1024 As tile_height
, tiles AS (
((ST_XMax(e) - ST_XMin(e))/dims.perpixel)::integer
, ((ST_YMax(e) - ST_Ymin(e))/dims.perpixel)::integer
, ST_XMin(e), ST_YMax(e)
, perpixel, -perpixel, 0, 0,f.srid
,tile_height/perpixel)::geometry As geom
FROM (SELECT ST_Extent(the_geom) As e, ST_SRID(the_geom) As srid
from neighborhoods GROUP BY srid) As f
CROSS JOIN dims)
SELECT 'gdal_translate -of jpeg big_o_raster.sid -projwin '
|| ST_XMin(tiles.geom) || ' ' || ST_YMax(tiles.geom)
|| ' ' || ST_XMax(tiles.geom) || ' ' || ST_YMin(tiles.geom)
|| ' tile_' || lpad(ROW_NUMBER() OVER()::text,5, '0') || '.jpg' As gdal_command
FROM tiles INNER JOIN
(SELECT ST_Union(the_geom) AS geom FROM neighborhoods) As f
ON ST_Intersects(tiles.geom, f.geom);
The output of the above sql are gdal_translate commands of the form:
gdal_translate -of jpeg big_o_raster.sid -projwin 772291.910466813 2970042.61545891 773311.910466813 2969022.61545891 tile_00001.jpg
gdal_translate -of jpeg big_o_raster.sid -projwin 786571.910466813 2970042.61545891 787591.910466813 2969022.61545891 tile_00002.jpg
gdal_translate -of jpeg big_o_raster.sid -projwin 787591.910466813 2970042.61545891 788611.910466813 2969022.61545891 tile_00003.jpg
gdal_translate -of jpeg big_o_raster.sid -projwin 788611.910466813 2970042.61545891 789631.910466813 2969022.61545891 tile_00004.jpg