[GIS] Clipping raster by shapefile in PyQGIS

clippyqgis

Currently I am working with some satellite images (grib files) in pyqgis. As I want to work only with the cells in my study area, I need to clip the raster by a shapefile. I know I can do this by Raster tools in QGIS, but since I have a lot of rasters to work with, I wanted to use python and save time.

Perhaps somebody knows the pyqgis commands to do this?

I am new to QGIS and python.

Best Answer

The raster tool of QGIS uses GDAL commands

You can also use these GDAL commands (gdalwrap or gdal_translate) in a Python script, look at how to for loop a folder to batch clip rasters by polygon using python and QGIS or Extract raster windows from satellite imagery with GDAL and in pure Python Python GDAL/OGR Cookbook: Python Clip a GeoTiff with Shapefile

But you can do the same thing without PyQGIS or GDAL with the Python module PyShp: (Shapefile): GeospatialPython: Clip a Raster using a Shapefile

And you can also use Fiona and rasterio -> Warping images with rasterio

Related Solutions

[GIS] Clipping raster brick object in R

In order to write a raster stack or brick, the Interleave band can be used in the option argument.

    writeRaster(ndvi_stack, filename="<filename>.tif", options="INTERLEAVE=BAND", overwrite=TRUE)

[GIS] pyQGIS Merge raster layers

You don't need processing to execute gdal_tools you can use subprocess module to execute your command. Even Processing use the subprocess so why don't you ? here a function to use it :

import sys
from subprocess import Popen, PIPE, STDOUT
# Create a class herit from Exception class
class ProcessException(Exception):
    pass
# function to run the command
def execute(command):

    process = Popen(command.split(" "), shell=True, stdout=PIPE, stderr=STDOUT, stdin=PIPE)
    # Poll process for new output until finished
    while True:
        nextline = process.stdout.readline()
        if nextline == '' and process.poll() is not None:
            break
        # You can add here some progress infos

    output = process.communicate()[0]
    exitCode = process.returncode

    if (exitCode == 0):
        return output
    else:
        raise ProcessException(command, exitCode, output)

To use it :

try:
    execute(u"gdal_merge.py -n 0 -a_nodata 0 -of GTiff -o /path/of/your/output/mozaik.tif /path/of/your/input/terrain_corrected01.tif /path/of/your/input/terrain_corrected02.tif")
except ProcessException as p:
    print(u"something goes wrong : {}".format(p))

--EDIT-- Where is those command

Gdal comes with QGIS as you install QGIS. In Windows they are in the :

C:\path\of\OSGEO4W\bin\

In linux it's in the :

/usr/bin/
# for the manual
/usr/share/man/manl/gdal_merge.l.gz

In my opinion if you don't need the QGIS API it will be faster and easy to use it in a bash script. But it depends on your need.

Best Answer

Related Solutions

[GIS] Clipping raster brick object in R

[GIS] pyQGIS Merge raster layers

Related Question