Clip raster by shapefile in parts?

Question

I am trying to write .sh script command to clip raster by polygonal shapefile. Shapefile have name field with unique name of polygon feature. On the output must be rasters with name of polygons in shapefile. I'm use Linux, GDAL. Now I use command

gdalwarp -dstnodata -9999 -cutline data/shapefile.shp data/input.tiff data/output.tiff

but in output I receive one raster.

Example on the picture:

Answer 1

Reorganise your shapefile so that one shapefile contains one feature (A,B,C in your case) only

Then use a loop like

for i in A B C; do
  gdalwarp -cutline $i.shp ... $i.tif
done

to create each output raster.

Example of script:

#!/bin/sh
#    "shp" - folder for shapefiles 
#    "outputraster" - folder for output rasters
cd /home/user/cliprasters/
#  delete previous files
rm ./outputraster/*.tiff
#  run
for fname1 in ./shp/*.shp do
  fname2="$(basename $fname1 .shp)"
  gdalwarp -dstnodata -9999 -cutline ./shp/$fname2.shp -crop_to_cutline -of GTiff ./input.tiff ./outputraster/$fname2.tiff
done

Answer 2

A simpler solution, although less elegant is simply use rasterio, gdal and numpy. Let's say that in our shapefile we have a column called 'column' containing the features names (our polygons), let's say that our polygons are 'A', 'B' and 'C'.

import gdal, os, rasterio
import numpy as np

im = 'image.tif'
shp = 'shapefile.shp'

#open the image using rasterio and reshape it to a np.array
a = rasterio.open(im).read()
a = a.reshape((a.shape[1],a.shape[2]))

# gdal is weird when operating with folders sometimes, 
#    so we just move to the folder where the image is
os.chdir(folder) #folder with image

#now we run gdalwarp with the appropriate parameters in a loop 
#  over our polygons. This is the part that could be improved
# by using fiona or geopandas, something to extract the features
# names and dump them in a list.
# The trick with the features names is to use single AND double quotes, like "'this'", because gdal wants to get something like 'this' as a parameter.


polnames = ["'A'","'B'","'C'"] 
os.chdir(folder)
for p in polnames:
    os.system('gdalwarp -cutline %s -cwhere "%s = %s" -crop_to_cutline %s %s_%s' %(shp,c,p,im,p,im)) 

Now you have the image cut in pieces according to your features in the shapefile. Then, we just open the images and get the data.

# open one image using rasterio
b = rasterio.open('Aim.tif').read().reshape((b.shape[1],b.shape[2]))
# if you want, you can flatten and have a 1d array
b = b.flatten()

Done! I said it is not elegant, there is a lot of room to improve, but it works.