16

I'm writing a simple utility to crop batches of multi-band geotiff raster files to the same (smaller) area. Using gdalwarp, I can easily crop a file using a single-polygon clipping shapefile:

gdalwarp -cutline clipper.shp -crop_to_cutline input.tif output.tif

However, the actual area I want to clip to will always be initially defined by another geotiff raster file, not a shapefile. It would be nice if I could use the extent of that raster as the clipping file, but I'm not sure how to do this. Unsurprisingly, the following doesn't work (it doesn't raise an error, it just doesn't produce anything):

gdalwarp -cutline clipper.tif-crop_to_cutline input.tif output.tif

So, my question is, is there a way to supply a raster to gdalwarp -cutline? Alternately, is there another gdal function that can clip a raster using another raster? If neither of these are possible, is there a very simple way to produce a shapefile with a single polygon defined by the extent of a raster?

This code will be wrapped in a more extensive python script, so I can use either command line gdal utilities or any of the python bindings for gdal.

As a side note, I know that I could easily just make a clipping shapefile that covers the extent of my raster in QGIS. I may wind up doing that if I don't find a straightforward solution, but I will ultimately wind up using this utility on dozens if not hundreds of areas as part of a large automated analysis, so I'd prefer not to have a tedious manual step even if it is very easy.

2
  • This solution only works if they are in the same SRS. Works fine for me. U Can try put "MEM" in the output and list bands from a stack.
    – Fronza
    Apr 7, 2020 at 14:51
  • @Fronza , hey can you put the request that works for you , thanks in advance
    – work
    Jun 16, 2021 at 10:16

4 Answers 4

12

I don't know if it's possible to clip a raster with an other raster but you could use gdaltindex to build the shapefile with the extent of your raster.

http://www.gdal.org/gdaltindex.html

2
  • 4
    gdaltindex works perfectly to make a clipping shapefile from my initial raster. To solve the problem I use gdaltindex clipper.shp clipper.tif, followed by gdalwarp -cutline clipper.shp -crop_to_cutline input.tif output.tif
    – Joe
    Dec 12, 2014 at 17:03
  • 1
    I was using this approach but found that it sometimes was a single pixel off in the clipped version. I think it's more straightforward to compute your target extents a la Xavier's answer below and then use gdalwarp and specify -te_srs to handle mismatched CRSs.
    – Jon
    Dec 18, 2019 at 5:22
9

For irregular polygons, and assuming that your geotiff raster file is a binary raster, you could use GDAL_Calc:

GDAL_Calc.py -A Mask.tif -B CutBigImageToClip.tif --outfile=SmallerFile.tif --NoDataValue=0 --Calc="B*(A>0)" 

This query will populate 0 where Mask.tif <= 0 and BigImage where the Mask > 0. To do this both rasters must be the same cell size, rows and columns. To extract the same extents use GDAL_Translate with the -projwin ulx uly lrx lry option (box is in projected coordinates), but ensure that the projwin box does not extend over the edges of either raster.

GDAL_Translate -of GTIFF -projwin ulx uly lrx lry BigImageToClip.tif CutBigImageToClip.tif

Substitute values for the projwin box derived from the Mask.

1
  • 1
    +1 This is useful information, but I think I can solve my problem in fewer steps using @lejedi's answer.
    – Joe
    Dec 12, 2014 at 17:00
5

The solution in Python directly, without making shape:

import gdal
from gdalconst import GA_ReadOnly

data = gdal.Open('img_mask.tif', GA_ReadOnly)
geoTransform = data.GetGeoTransform()
minx = geoTransform[0]
maxy = geoTransform[3]
maxx = minx + geoTransform[1] * data.RasterXSize
miny = maxy + geoTransform[5] * data.RasterYSize
call('gdal_translate -projwin ' + ' '.join([str(x) for x in [minx, maxy, maxx, miny]]) + ' -of GTiff img_orig.tif img_out.tif', shell=True)
2
  • 1
    NB: This solution only works if they are in the same SRS.
    – Skylion
    Nov 16, 2016 at 18:10
  • 1
    @Skylion But you can easily account for this by including the -te_srs option, although you also need to gdalwarp instead with the -te option.
    – Jon
    Dec 18, 2019 at 5:20
0

Here u can use this code to run an intersection raster tool. You need to be aware of CRS and pixel size (cols and rows).

# This file is part of Brazil Data Cube Validation Tools.
# Copyright (C) 2020.

# Python Native
import os
# 3rd party
import gdal
import numpy


def raster_intersection(ds1, ds2, nodata1=None, nodata2=None, output_name1=None, 
output_name2=None):
"""Perform image intersection of two rasters with different extent and 
projection.
    Args:
        ds1 (GDAL dataset) - GDAL dataset of an image
        ds2 (GDAL dataset) - GDAL dataset of an image
        nodata1 (number) - nodata value of image 1
        nodata2 (number) - nodata value of image 2
        output_name1 (string) - path to output intersection of ds1
        output_name2 (string) - path to output intersection of ds2
    Returns:
        dataset1 (GDAL dataset), dataset2 (GDAL dataset): intersection dataset1 
and intersection dataset2.
"""
###Setting nodata
nodata = 0
###Check if images NoData is set
if nodata2 is not None:
    nodata = nodata2
    ds2.GetRasterBand(1).SetNoDataValue(nodata)
else:
    if ds2.GetRasterBand(1).GetNoDataValue() is None:
        ds2.GetRasterBand(1).SetNoDataValue(nodata)

if nodata1 is not None:
    nodata = nodata1
    ds1.GetRasterBand(1).SetNoDataValue(nodata1)
else:
    if ds1.GetRasterBand(1).GetNoDataValue() is None:
        ds1.GetRasterBand(1).SetNoDataValue(nodata)

### Get extent from ds1
projection = ds1.GetProjectionRef()
geoTransform = ds1.GetGeoTransform()

###Get minx and max y
minx = geoTransform[0]
maxy = geoTransform[3]

###Raster dimensions
xsize = ds1.RasterXSize
ysize = ds1.RasterYSize

maxx = minx + geoTransform[1] * xsize
miny = maxy + geoTransform[5] * ysize

###Warp to same spatial resolution
gdaloptions = {'format': 'MEM', 'xRes': geoTransform[1], 'yRes': 
geoTransform[5], 'dstSRS': projection}
ds2w = gdal.Warp('', ds2, **gdaloptions)
ds2 = None

###Translate to same projection
ds2c = gdal.Translate('', ds2w, format='MEM', projWin=[minx, maxy, maxx, miny], 
outputSRS=projection)
ds2w = None
ds1c = gdal.Translate('', ds1, format='MEM', projWin=[minx, maxy, maxx, miny], 
outputSRS=projection)
ds1 = None

###Check if will create file on disk
if output_name1 is not None or output_name2 is not None:
    driver = gdal.GetDriverByName("GTiff")
    if output_name1 is None:
        output_name1 = 'intersection1.tif'
    if output_name2 is None:
        output_name2 = 'intersection2.tif'
else:
    driver = gdal.GetDriverByName("MEM")
    output_name1 = ''
    output_name2 = ''

dataset1 = driver.Create(output_name1, xsize, ysize, 1, 
ds1c.GetRasterBand(1).DataType)
dataset1.SetGeoTransform(geoTransform)
dataset1.SetProjection(projection)
dataset1.GetRasterBand(1).SetNoDataValue(nodata) ###Setting nodata value
dataset1.GetRasterBand(1).WriteArray(ds1c.GetRasterBand(1).ReadAsArray())

dataset2 = driver.Create(output_name2, xsize, ysize, 1, 
ds2c.GetRasterBand(1).DataType)
dataset2.SetGeoTransform(geoTransform)
dataset2.SetProjection(projection)
dataset2.GetRasterBand(1).SetNoDataValue(nodata) ###Setting nodata value
dataset2.GetRasterBand(1).WriteArray(ds2c.GetRasterBand(1).ReadAsArray())

ds1c = None
ds2c = None

return dataset1, dataset2


def raster_absolute_diff(ds1, ds2, nodata1=None, nodata2=None, 
output_file=None):
"""Perform image absolute difference (support different extent and projection).
    Args:
        path1 (string) - path to image 1 (reference)
        path2 (string) - path to image 2 (target)
        output_dir (string) - path to output files
        nodata1 (number) - nodata value of image 1
        nodata2 (number) - nodata value of image 2
    Returns:
        dataset (GDAL dataset): dataset containing absolute difference between 
ds1 and ds2.
"""
if output_file is None:
    output_file = 'abs_diff.tif'
ds1_intersec, ds2_intersec = raster_intersection(ds1, ds2, nodata1, nodata2, 
None, None)

### Read bands with numpy to algebra
nodata = ds1_intersec.GetRasterBand(1).GetNoDataValue()
bandtar = numpy.array(ds1_intersec.GetRasterBand(1).ReadAsArray().astype(float))
fill_bandtar = numpy.where(bandtar == nodata)
bandref = numpy.array(ds2_intersec.GetRasterBand(1).ReadAsArray().astype(float))
fill_bandref = numpy.where(bandref == nodata)

### Get extent from ds1
projection = ds1.GetProjectionRef()
geoTransform = ds1.GetGeoTransform()
[cols, rows] = ds1.GetRasterBand(1).ReadAsArray().shape

ds1 = None
ds2 = None
diff = numpy.abs(bandtar - bandref)
diff[fill_bandtar] = nodata
diff[fill_bandref] = nodata

###Check if will create file on disk
if output_file is not None:
    driver = gdal.GetDriverByName("GTiff")
else:
    driver = gdal.GetDriverByName("MEM")
    output_file = ''

dataset = driver.Create(output_file, rows, cols, 1, 
ds1_intersec.GetRasterBand(1).DataType)
dataset.SetGeoTransform(geoTransform)
dataset.SetProjection(projection)
dataset.GetRasterBand(1).SetNoDataValue(nodata)
dataset.GetRasterBand(1).WriteArray(diff)

return dataset
1
  • 2
    Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.
    – Community Bot
    Mar 22 at 19:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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