3

I'm looking for a way to clip a tif with a shapefile containing multiple polygons like here.

When looking around all the options I find are really long like here and many similar ones. However I can't believe that there is no one line solution for this. Am I missing something?

4
  • 1
    Do you need a pure python solution?
    – YoLecomte
    Oct 17, 2017 at 15:20
  • No, but I am limited in what I can and cannot install on the machine I'm using
    – ArnJac
    Oct 17, 2017 at 15:31
  • You can try my answer which use gdal
    – YoLecomte
    Oct 17, 2017 at 15:46
  • I agree that it is kind of a hassle if you want to stay at the gdal level. You can use rasterio, though, as shown here: automating-gis-processes.github.io/CSC18/lessons/L6/… , but you'd need to alter that to get the bounding box coordinates of your shapefile, which is very easy.
    – Jon
    May 19, 2018 at 21:07

2 Answers 2

1

if gdal is possible you can use gdalwrap with the -cutline option (see http://www.gdal.org/gdalwarp.html)

You can call a gdal command in python with subprocess (example here : http://chris35wills.github.io/subprocess_gdal/)

but I prefer create the command line as a list (to keep control on the option and play with variable) and then join it like that :

cmd =  " ".join(['gdal_translate','-co','"COMPRESS=LZW"', input_folder+file,output_folder+file])
3
  • I still need to install gdal to do this right? And this is not the one that I can install through pip I assume? The machine doesn't let me install through the exe of osgeo4...
    – ArnJac
    Oct 17, 2017 at 16:01
  • I'm not a specialist of installation and requirement and this is not the subject of this question. You could ask another question for that specific problem "installing gdal properly on a your OS machine...". Maybe, the solution through pip could work!
    – YoLecomte
    Oct 17, 2017 at 16:05
  • I recommend Anaconda to manage environments and package dependencies. Installing gdal would be conda install gdal. I would think that the pip install gdal would be fine also.
    – Jon
    May 19, 2018 at 21:09
1

Here's a script I use to clip rasters by shapefiles--note that it returns a raster object that contains the array, masking, and georeferencing. It would be another step to write the array as a geotiff, but not a difficult one. Note that you need to install the fiona and rasterstats packages.

import fiona
from rasterstats.io import Raster

def clip_raster_by_shp(rasterpath, shppath, band=1):
    """
    Given an input raster (path) and shapefile (path), returns the raster object
    corresponding to the raster defined by the bounds of the shapefile
    """

    with fiona.open(shppath) as f:
        geom_bounds = f.bounds

    with Raster(rasterpath, band=band) as raster_obj:
        rast = raster_obj.read(bounds=geom_bounds) 

    return rast

To see what the rast object contains, after running this type rast. then press tab (this assumes you're using an IDE like Spyder). But if you're not manipulating the clipped raster and you need to save it as a geotiff, I would definitely go with the other answer.

1
  • Thanks so much, this really helped me!
    – edge-case
    Sep 2, 2019 at 14:16

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.