6

I am trying to split a single-band raster (land cover) into smaller parts representing administrative units which are into another shapefile, containing their name as an attribute.

U I did this using SAGA "clip raster with polygon" (in QGIS 2.18) to split the raster and then naming the resulting rasters one by one (by hand) with the name of the administrative unit.

enter image description here

The problem is, I now need to do this process over ±2000 polygons, so doing it by hand is the last resort options.

I tried different options:
a) Using SAGA "clip raster with polygon", iterate over the polygon layer. Works good, but cannot specify an attribute to use a name.

b) Batch-process SAGA "clip raster with polygon" using individual administrative polygons, named after their administrative ID. Works fine, the name of the resulting rasters is as I want it, but the process is long and boring (batch-processing with 2000 rows, filled in part by hand, see image).

enter image description here

So, is there a way of speeding up the process? Is it possible to iterate over a layer, specifying an attribute to use as a name for the resulting layers? Is it possible to auto-fill the batch-processing window with layers open in the "layers panel"?

2
  • Possible duplicate of Clip raster by shapefile in parts?
    – ImanolUr
    Commented Aug 27, 2018 at 15:21
  • Have you an automaticly method to add raster to each polygone, like adding polygones by select all. because Enter the name of the raster by copy-pasting it on every row take all of time to import it to each row?
    – resuk
    Commented Jan 14, 2022 at 16:44

2 Answers 2

6

Your question is quite similar to this one, but it does not say how to automate the separation of the files, so I will answer:

import glob
import fiona
from subprocess import Popen

with fiona.open('polygons.shp', 'r') as dst_in:
    for index, feature in enumerate(dst_in):
        with fiona.open('separated/polygon{}.shp'.format(index), 'w', **dst_in.meta) as dst_out:
            dst_out.write(feature)

polygons = glob.glob('separated/*.shp')  ## Retrieve all the .shp files

for polygon in polygons:
    feat = fiona.open(polygon, 'r')
    name = feat['properties']['Name_of_your_attribute']  
    command = 'gdalwarp -dstnodata -9999 -cutline {} ' \
              '-crop_to_cutline -of GTiff ./input.tiff ./outputraster/{}.tiff'.format(polygon, name )
    Popen(command, shell=True)

This code separates the polygons to different files and then clips the raster with GDAL

6

Hello and thanks for your code ImanoIUr.

As I am not familiar (at all) with python, I tried to find a not-python way of doing this task. I came up with this, and it works fine, except it is more time-consuming than automating it with python. I am posting this answer because it could be useful to other people that do not use python.

  1. Split the polygon (administrative) layer with "split polygon layer", using the attribute corresponding to the name you want to give your end-result rasters.
  2. Batch-process SAGA "clip raster with polygon". Enter the name of the raster by copy-pasting it on every row (image below). For the "polygons" column, click the (...) and "select by filesystem", then go to your splitted shapefile directory and select all your shapefiles at once (image below). For the "clipped" column, click (...), specify the folder you want the output rasters to go, an click OK. You will be prompted with the "autofill settings" box, for wich you select "polygons" as the naming parameter (see other image).
    enter image description here
    enter image description here
  3. Hit "run" and voilà! The output folder is filled with rasters clipped with each polygon of your original polygon layer. Assuming you named them right (in step 1), the name of the administrative unit will appear as a part of your end result raster's name.

Although not perfect, as the name of the final product is not exactly what I wanted it to be, the info I needed is there; which make this solution good for me.

Your Answer

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

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