Clipping GeoTIFF with shapefile?

Question

I have a GeoTIFF file and a shapefile. I want to clip the GeoTIFF file using the polygons defined in the shapefile in such a way that it produces a separate output for each polygon (which is different from the gdalwarp usecase, which produces a single output by default).

Furthermore, I want to do this in an efficient manner. Iterating over each polygon and calling gdalwarp will reload the GeoTIFF each and every time, which is something I want to avoid.

If some tool (or some options of gdalwarp) exists to solve this use case, it would be very interesting to know.

EDIT: The piece of code that we are currently using:

def cut_region(fname, ids, output_folder, input_pic, input_shape):
    photo_name  = input_pic + fname + ".jp2"
    shape_name  = input_shape + fname + ".shp"
    raster = gdal.Open(photo_name)
    projection = raster.GetProjectionRef()

    for id in ids:
        output_name = output_folder + fname + '/' + str(id) + ".tif"
        output = gdal.Warp(output_name,
           raster,
           format='GTiff',
           cutlineDSName=shape_name,
           cutlineSQL="SELECT * from " + fname + " where ID='" + str(id) + "'",
           dstSRS = projection,
           cropToCutline=True,
           multithread=True,
           dstAlpha=True
        )

Answer 1

In QGIS, you can use the iterator of the clip tool to iterate over the polygon shapefile, and it will clip the GeoTiff raster based on each polygon inside the polygon shapefile.

For example, the following image shows a DEM raster and a polygon shapefile overlaying the raster file. This shapefile is composed of several grid polygons:

Using Clip raster by mask layer, which you can find it from Processing toolbox -> GDAL -> Vector geoprocessing -> Clip raster by mask layer, you can use the iterator (the green arrow) beside the polygon shapefile:

It will produce a separate clipped image based on each polygon inside the polygon shapefile:

Please note that you need to select the data type to be similar to data type of the original GeoTiff raster.

I used QGIS 3.2.2, but it is also available in QGIS 2.18.23.

Answer 2

I have found an efficient solution using a combination of rasterio and rasterstats. At this point it removes all georeferences (which is not required for our purposes), but these shouldn't be hard to add back in.

import fiona
import rasterio
import shapely.geometry
import numpy as np

from rasterstats.io import Raster
from PIL import Image


tif_filename = 'sample_sat.tif'
with Raster(tif_filename, band=1) as raster_obj_1:
    with Raster(tif_filename, band=2) as raster_obj_2:
        with Raster(tif_filename, band=3) as raster_obj_3:
            index = 0
            for feat in fiona.open('sample.shp'):
                polygon_geometry = feat['geometry']
                polygon = shapely.geometry.Polygon(polygon_geometry['coordinates'][0])
                polygon_bounds = polygon.bounds

                raster_subset_1 = raster_obj_1.read(bounds=polygon_bounds)
                polygon_mask = rasterio.features.geometry_mask(geometries=[polygon_geometry],
                                                    out_shape=(raster_subset_1.shape[0],raster_subset_1.shape[1]),
                                                    transform=raster_subset_1.affine,
                                                    all_touched=False,
                                                    invert=True)

                raster_subset_2 = raster_obj_2.read(bounds=polygon_bounds)
                raster_subset_3 = raster_obj_3.read(bounds=polygon_bounds)

                masked_1 = raster_subset_1.array * polygon_mask
                masked_2 = raster_subset_2.array * polygon_mask
                masked_3 = raster_subset_3.array * polygon_mask

                masked_all = np.dstack([masked_1, masked_2, masked_3])

                img = Image.fromarray(masked_all[:, :, :].astype('uint8'), 'RGB')
                img.save('out/' + str(index) + '.jpg')
                index += 1

Answer 3

The following example demonstrates how you can replicate this functionality from the command-line (which might be even more efficient).

Let us assume that you have a set of polygons (e.g. poly_extents.shp), which is located within your raster (input_image.tif). Each polygon would also have a unique field called id numbered from 1:5 (see image below).

You could then use a bash for loop to iterate over each unique polygon, which would be passed to gdalwarp using the SQL command defined using the -csql option.

for id in {1..5}
     do echo $id  
     gdalwarp  -cutline poly_extents.shp -crop_to_cutline \ 
               -csql "SELECT * FROM poly_extents WHERE id = $id"  \ 
               input_image.tif output_clip_${id}.tif
done

The result would be a set of 5 GeoTiffs: