Clip By Polygon, Clip By Mask, Clip By File

[maphew]

A feature I wish for almost every week is the ability to clip an image by an irregular polygon. I envision it working something like this:

1. Clip By Polygon
2. Clip by Raster
3. Clip By File
4. Globals

Clip By Polygon

gdal -clip poly.shp infile.tif outfile.tif  {options}

Cookie cutter an image. Outside the polygon is nodata, inside is input.tif values. Output extent is set to shapefile projwin extents.

gdal -clip poly.shp -reverse in.tif out.tif

Punches a hole in the image. In output raster interior of polygon is nodata. Output extent is the larger of raster or shapefile.

gdal -clip poly.shp -multi -name [attrib] in.tif outdir

Create tiles. For each polygon in poly.shp create a raster clipped to extent of that polygon. Use [attrib] field in poly.shp to name output rasters. If -name is omitted, sequentially number output (in.tif --> outdir/in_00.tif).

General

If clipping shapefile consists of multiple polygons, clipping boundary becomes outer boundary of all polys (merge polygons before clip). Polgons can be discontiguous (voids and islands are respected).

Clip by Raster

gdal -clip mask.tif -mask 255,0,0 in.tif out.tif

Same as clip by poly, but treat RGB value of 255,0,0 as the nodata area (e.g. is the outside of the clipping polygons). Output extent is shrunk/expanded to range of data areas.

gdal -clip mask.tif -mask 255,0,0 -reverse in.tif out.tif

Punch a hole. Same as clip by poly reversed. Output extent is larger of inputs.

General

If no mask is specified, clipping boundary is set to extents of clipping image. Nodata in clipping image is nodata in output image.

Clip By File

gdal -clip points.shp -extent-only in.tif out.tif
gdal -clip lines.shp -extent-only in.tif out.tif
gdal -clip polys.shp -extent-only in.tif out.tif
gdal -clip input.jp2 -extent-only in.tif out.tif

Ignore all boundary/mask logic and clip raster to projwin extent of any ogr/gdal supported data source.

Globals

Standard gdalwarp and gdal_translate options are respected, e.g. -co compress=lzw, -of jp2kak, -t_srs +utm=11, -outsize 50% 50%, etc.

[maphew]

whups, the internal table of contents is wrong. It's supposed to refer to each of the subheadings below. I can't remove it though, sorry.

[laurent]

Do this feature have a chance to be released ? Maybe it's just a dream ?

[warmerdam]

Laurent,

This is currently a request enhancement but no one has indicated a willingness to implement it. So effectively, yes, it is a dream.

Note that gdal_rasterize can already do some of this (especially with the new -i flag - similar to -reverse above). But gdal_rasterize requires an output file to pre-exist. I think the best approach to this might be further work on gdal_rasterize.

[maphew]

...I implemented part of what is asked in that ticket in the just released Mirone 1.3.0 We can now clip an image by an arbitrary polygon and create masks from the ensemble of plotted polygons (just right click on polygon and select what you want).

Joaquim Luis -- ​http://www.nabble.com/updates-on-masking-raster-data-tp17028508p17028831.html

Mirone's homepage: ​http://w3.ualg.pt/~jluis/mirone/

added to this ticket on the chance there may be an opportunity to share code between the projects for this feature.

[Mike Taves]

It would be very convenient to have a "-clip" option available for gdalwarp and gdal_translate (I'm not sure the scope of this enhancement, but this is the closest ticket that I could find to add to).

I'd like to also suggest a "Clip by WKT" feature to this enhancement, for example:

-clip "POLYGON((524813.7 6522424.5,524829.9 6522424.5,524829.9 6522400.3,524813.7 6522400.3,524813.7 6522424.5))"

which would use the same projection as the input raster.

Or, by specifying a different projection for the mask shape, such as Lat/Long WGS84:

-clip "SRID=4326;POLYGON((-126.5 59.25,-126.5 52.5,-126.25 52.5,-126.25 59.25,-126.5 59.25))"

[bicealyh]

I want to know does GDAL support clip an image by an irregular polygon,now?

[warmerdam]

The situation is unchanged. This feature has not been implemented, but -srcwin (or -projwin) in combination with gdal_rasterize can accomplish some of this in a rather manual way.

I prefer not to set unresourced enhancement requests above normal priority.

[Even Rouault]

Clip by WKT can be achieved by using the cutline option of gdalwarp and a small CSV file containing the WKT (trunk feature).

The following python script produces a white square.tif and a CSV file describing a polygon with a hole in the extent of square.tif

import gdal
import osr

ds = gdal.GetDriverByName('GTiff').Create('square.tif', 1000, 1000, 3)
ds.SetGeoTransform([500000, 1, 0, 4501000, 0, -1])
sr = osr.SpatialReference()
sr.ImportFromEPSG(32631)
ds.SetProjection(sr.ExportToWkt())
ds.GetRasterBand(1).Fill(255)
ds.GetRasterBand(2).Fill(255)
ds.GetRasterBand(3).Fill(255)
ds = None

f = open('small_square_with_hole.csv', 'wt')
f.write('WKT,dummy\n')
f.write('"POLYGON((500250 4500250,500750 4500250,500750 4500750,500250 4500750,500250 4500250),(500375 4500375,500625 4500375,500625 4500625,500375 4500625,500375 4500625))",\n')
f.close()

Then I do 'gdalwarp -wo "INIT_DEST=255,0,255" square.tif square_cut.tif -cutline small_square_with_hole.csv' and I get a nice extract of my white square with a hole and a background of pink.

[jason]

Has anyone ever got the cutline to work in C++?

I have it working off the command line no problem, but it doesn't seem to do anything when I populate that argument in C++ as below:

*Create a OGRGeometry variable OGRGeometry thisGeom; parse an input WKT to char char* wktIn = (char*) pIn_wkt_geom; Create and polulate thisGeom with the input WKT geometry OGRErr err = OGRGeometryFactory::createFromWkt(&wktIn, poSRS, &thisGeom);

Add thisGeom as the cutline argument psWarpOptions->hCutline = thisGeom;*

Has anyone achieved this or do you know if it is supported?

Cheers Jason

[Jukka Rahkonen]

Jason, hardly anybody reads 8 years old tickets like this. Ask your question in gdal-dev mailing list.

[Even Rouault]

Removing obsolete milestone

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.