Best way to clip a large raster ECW file ?

Question

I am trying to clip a large ECW (details below) but the raster file is too big to be processed entirely.

Some details of the ECW below

Driver: ECW/ERDAS Compressed Wavelets (SDK 5.0)

file size: 50gb Size is 450000, 565081 Pixel size: 0.15 0.15 COLORSPACE=RGB COMPRESSION_RATE_TARGET=9 VERSION=2 Number of bands: 4

The area I want to clip is roughly a 1/5 of the original file.

Here are the methods I have tried out with no success:

1. Used Arcgis to save the ecw to tiff / other formats... (I quickly gave up)

2. Used Qgis and its clipper tool... the file creation stayed stuck at about 40%.

3. Used gdal_translate out of OSGeo4W with other options than Qgis. (Tried that thinking that MAYBE freeing some memory not using Qgis would the trick)

4. Used gdal_retile thinking that I would cut the image into pieces and grab the one I wanted. The command "gdal_retile -ps 10000 10000 -of ecw -tileIndex tile.shp -targetDir input.ecw This crashed even quicker"

Does anyone have an idea?

For info I run windows 7 64bits on an i5-3470 3.2Ghz with 16gb of ram.

Answer 1

Other ideas you could try:

1. gdal_translate with the -srcwin switch
2. gdalwarp with -cutline and -crop_to_cutline and -wm switches. The last one specifies memory for caching and may get you over the issues you had using clipper in QGIS (as this is essentially the same function)
3. QGIS raster calculator setting the extent to the area you want (simpler than clipper).
4. SAGA->Clip grid with polygon - who knows, it might be more memory efficient.
5. Code a solution using Python and Numpy/SciPy to read just a subset of the raster into memory and save it.

I suspect that doing a very simple crop to extent (either with the raster calculator or -srcwin switch in gdal_translate) will be less memory hungry than cropping with a polygon because you have no geometry checks and conversions. Option 5 should use the least memory as you are only reading in what you need. Have a look at this tutorial if you need a 'how to' (adapt the bit on reading and writing by block).

Answer 2

You can cut it directly with gdal's tool gdal_translate if you know the coordinates of your Area of Interest, if its georeferenced:

gdal_translate -projwin [ulx uly lrx lry] infile outfile  

If not use the -srswin flag like this:

gdal_translate -srcwin [xoff yoff xsize ysize] infile outfile.

Another option is to build a 'virtual' raster (of a few kilobytes) that point at your initial dataset, with gdalbuildvrt.

gdalbuildvrt -te [xmin ymin xmax ymax] infile outfile.vrt

Unfortunally you will need a georeferenced file for using gdalbuildvrt.

After you clip, don't forget to create pyramids for easy viewing. The following command will build external pyramids compressed with the DEFLATE (lossless) algorithm:

gdaladdo -ro --config COMPRESS_OVERVIEW DEFLATE outfile 2 4 8 16

As a final step your can calculate the statistics as well to avoid some stupid problems with a particular commercial program:

gdalinfo -stats outfile

Answer 3

A straight 'clip' using gdalwarp should work (I know this is a hella-old question: 18 months IRL is like a geological epoch in internet years).

I have a 70Gb aerial (ECW, 94000x81000 pixel at 10cm/px), and GDAL can arbitrarily clip it with a shapefile using

gdalwarp -cutline [clipfile] -crop_to_cutline [infile] [outfile]

at the Windows command-line. (I realise that this solution requires building a shapefile for the target area of interest, but that's not the greatest challenge known to GIS).

To extract a half-suburb-sized chunk takes ~0.4sec; to slice the file into quarters takes 4sec. My machine's specs are not dramatically different (i7-4770 @ 3.4GHz, 16GB RAM, Win7-64 Ultimate).