[Gdal-dev] Masking portions of a raster outside of a vector boundary

[Frank Warmerdam]

Avi A Blackmore wrote:
> Hi,
> 
> 	I'm trying to determine the best means of "masking out" the parts of a 
> raster image that lie outside of a vector boundary, setting those areas to 
> a null value.  Are there any existing GDAL or OGR tools which could be 
> used to do this?
> 
> 	If not, any suggestions for a general algorithm using GDAL/OGR?

Avi,

The gdal_rasterize utility can mask out areas *inside* polygons.  If you can
wrap a "universe polygon" around all your polygons then that would give you
what you want by inverting inside and outside.

I have wanted to be able to do this "mask outside" sort of operation for
a long time, so I have finally added the ability to gdal_rasterize.  It
now (in subversion as of a few minutes ago) has a -i (invert) flag that
will do this for you.  So for instance:

   gdal_rasterize -burn 128 utm.shp -i -l utm out.tif

will mark all areas outside the polygons in utm.shp to a value
of 128 in out.tif.

Keep in mind that "per feature" burn values no longer make sense in this
context.  The code will now will use the burn value associated with the
very first feature when in -i mode.

For those interested, the inverter merges all the geometries into one
geometry collection and then add a universe polygon as shown below.
A use of the much hated OGR C geometry api, and one of the few places
where I use std::vector.

/************************************************************************/
/*                          InvertGeometries()                          */
/************************************************************************/

static void InvertGeometries( GDALDatasetH hDstDS,
                               std::vector<OGRGeometryH> &ahGeometries )

{
     OGRGeometryH hCollection =
         OGR_G_CreateGeometry( wkbGeometryCollection );

/* -------------------------------------------------------------------- */
/*      Create a ring that is a bit outside the raster dataset.         */
/* -------------------------------------------------------------------- */
     OGRGeometryH hUniversePoly, hUniverseRing;
     double adfGeoTransform[6];
     int brx = GDALGetRasterXSize( hDstDS ) + 2;
     int bry = GDALGetRasterYSize( hDstDS ) + 2;

     GDALGetGeoTransform( hDstDS, adfGeoTransform );

     hUniverseRing = OGR_G_CreateGeometry( wkbLinearRing );

     OGR_G_AddPoint_2D(
         hUniverseRing,
         adfGeoTransform[0] + -2*adfGeoTransform[1] + -2*adfGeoTransform[2],
         adfGeoTransform[3] + -2*adfGeoTransform[4] + -2*adfGeoTransform[5] );

     OGR_G_AddPoint_2D(
         hUniverseRing,
         adfGeoTransform[0] + brx*adfGeoTransform[1] + -2*adfGeoTransform[2],
         adfGeoTransform[3] + brx*adfGeoTransform[4] + -2*adfGeoTransform[5] );

     OGR_G_AddPoint_2D(
         hUniverseRing,
         adfGeoTransform[0] + brx*adfGeoTransform[1] + bry*adfGeoTransform[2],
         adfGeoTransform[3] + brx*adfGeoTransform[4] + bry*adfGeoTransform[5] );

     OGR_G_AddPoint_2D(
         hUniverseRing,
         adfGeoTransform[0] + -2*adfGeoTransform[1] + bry*adfGeoTransform[2],
         adfGeoTransform[3] + -2*adfGeoTransform[4] + bry*adfGeoTransform[5] );

     OGR_G_AddPoint_2D(
         hUniverseRing,
         adfGeoTransform[0] + -2*adfGeoTransform[1] + -2*adfGeoTransform[2],
         adfGeoTransform[3] + -2*adfGeoTransform[4] + -2*adfGeoTransform[5] );

     hUniversePoly = OGR_G_CreateGeometry( wkbPolygon );
     OGR_G_AddGeometryDirectly( hUniversePoly, hUniverseRing );

     OGR_G_AddGeometryDirectly( hCollection, hUniversePoly );

/* -------------------------------------------------------------------- */
/*      Add the rest of the geometries into our collection.             */
/* -------------------------------------------------------------------- */
     unsigned int iGeom;

     for( iGeom = 0; iGeom < ahGeometries.size(); iGeom++ )
         OGR_G_AddGeometryDirectly( hCollection, ahGeometries[iGeom] );

     ahGeometries.resize(1);
     ahGeometries[0] = hCollection;
}


-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org