Chopping rasters with gdal_translate

BostonGIS BostonGIS Home About Postgres OnLine Journal PostGIS in Action

Thursday, May 09. 2013

Posted by Regina Obe in gdal, gis, postgis_raster, spatialite, waiting_postgis_21

Comments (0)
Trackbacks (0)

Chopping rasters with gdal_translate

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.

In a nutshell, we're using PostGIS raster technology (ST_Tile function introduced in PostGIS 2.1) that we demonstrated in Waiting for PostGIS 2.1 ST_Tile to create a grid because PostGIS geometry doesn't have cool gridding function like SpatiaLite has :). SpatialLite tesselation. Perhaps in PostGIS 2.2 we'll see some of these SpatiaLite niceties. However ST_Tile does the trick fairly nicely and quickly. For this example took under 600 ms to generate 1524 rows of GDAL commands.

WITH dims As (
 SELECT 10 As perpixel -- 10x10units/pixel  prevent exceeding size limit of raster
    , 1024 AS tile_width --1024 units of srs (srid)
    , 1024 As tile_height --1024 units of srs (srid)
    )
  , tiles AS (
  -- cut the empty raster into tiles of tile_width, tile_height
  -- spatial units
   SELECT ST_Tile(
    -- make an empty raster that is same dimensions
     -- as extent of our geometry
      ST_MakeEmptyRaster(
      -- we divide so that we don't go beyond 65536 pixel width/height integer (limit of raster)
       ((ST_XMax(e) - ST_XMin(e))/dims.perpixel)::integer
        , ((ST_YMax(e) - ST_Ymin(e))/dims.perpixel)::integer
        , ST_XMin(e), ST_YMax(e)
        -- note we have -perpixel for y direction
        -- because raster 
        -- is upperleft and geometry is lowerleft
       , perpixel, -perpixel, 0, 0,f.srid 
       )
  ,tile_width/perpixel
  ,tile_height/perpixel)::geometry As geom  --convert raster to geometry
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
-- only generate geometry tiles where our original tiles intersect
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
:

Trackbacks

Trackback specific URI for this entry

No Trackbacks

Comments

Display comments as (Linear | Threaded)

No comments

Add Comment

Name
Email
Homepage
In reply to
Comment	E-Mail addresses will not be displayed and will only be used for E-Mail notifications. To prevent automated Bots from commentspamming, please enter the string you see in the image below in the appropriate input box. Your comment will only be submitted if the strings match. Please ensure that your browser supports and accepts cookies, or your comment cannot be verified correctly. Enter the string from the spam-prevention image above:
	Remember Information? Subscribe to this entry