2

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
        )
4
  • 1
    Are you sure that gdalwarp is inefficient? It does not need to read the whole image from the disk each time. Just take care that your source tiff is tiled, not striped one.
    – user30184
    Sep 26, 2018 at 11:51
  • Using the gdal Python library the raster image indeed only needs to be loaded once. But even then it takes a couple of seconds to clip a single polygon, which seems relatively slow to me, especially if you want to do this for a couple of million polygons. I know I can parallelize the operation in order to make this task more manageable, but I was just expecting it to perform much better than this. Do these timings seem plausible? Does it depend on the size of the input raster image? And are there any other contributing factors to this performance?
    – jelleve
    Sep 27, 2018 at 12:43
  • Those who are familiar with GDAL and Python might be able to help you with optimizing your code if you show what you do now.
    – user30184
    Sep 27, 2018 at 13:11
  • I have come to realize that my question was not entirely correct. Instead of geotiff files, we are actually using jpeg2000 files, which seem to be a lot slower to process.
    – jelleve
    Oct 3, 2018 at 11:51

3 Answers 3

6

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:

enter image description here

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:

enter image description here

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

enter image description here

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.

3

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
2
  • 1
    any improvement regarding the georeference ? Jan 12, 2021 at 8:01
  • 1
    It has been a long time since I've worked on this, but I would think that you can just use the CRS metadata from the original raster and update the extent (which is probably equal to the polygon's bounds).
    – jelleve
    Jan 13, 2021 at 9:59
2

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).

enter image description here

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: enter image description here

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.