I wrote the function which cutting geotiff image by bounding box.
First image is original. At second you can see result off my code. I don't use gdalwarp or any console utilities. But I have no idea how to cut by geojson file. Also I can use only GDAL and numpy modules.
Here is my code:
import os, sys
from osgeo import gdal,gdalconst,osr
def cut_by_bounding_box(min_x, max_x, min_y, max_y):
xValues = [min_x, max_x]
yValues = [min_y, max_y]
#
# Register Imagine driver and open file
#
driver = gdal.GetDriverByName('GTiff')
driver.Register()
dataset = gdal.Open(filename)
if dataset is None:
print 'Could not open ' + filename
sys.exit(1)
#
# Getting image dimensions
#
cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount
#
# Getting georeference info
#
transform = dataset.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]
#
# Computing Point1(i1,j1), Point2(i2,j2)
#
i1 = int((xValues[0] - xOrigin) / pixelWidth)
j1 = int((yOrigin - yValues[0]) / pixelHeight)
i2 = int((xValues[1] - xOrigin) / pixelWidth)
j2 = int((yOrigin - yValues[1]) / pixelHeight)
new_cols = i2 - i1 + 1
new_rows = j2 - j1 + 1
#
# Create list to store band data in
#
band_list = []
#
# Read in bands and store all the data in bandList
#
for i in range(bands):
band = dataset.GetRasterBand(i+1) # 1-based index
data = band.ReadAsArray(i1, j1, new_cols, new_rows)
band_list.append(data)
new_x = xOrigin + i1 * pixelWidth
new_y = yOrigin - j1 * pixelHeight
new_transform = (new_x, transform[1], transform[2], new_y, transform[4], transform[5])
#
# Create gtif file
#
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, new_cols, new_rows, 3, gdal.GDT_Byte)
#
# Writting output raster
#
for j in range(bands):
data = band_list[j]
dst_ds.GetRasterBand(j+1).WriteArray(data)
#
# Setting extension of output raster
#
dst_ds.SetGeoTransform(new_transform)
wkt = dataset.GetProjection()
#
# Setting spatial reference of output raster
#
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#
# Close output raster dataset
#
dataset = None
dst_ds = None
if __name__ == '__main__':
# Imput/output file name and set directory
os.chdir('/home/sant/test/satellite_images')
filename = '20160501.tif'
output_file = '/home/sant/test/20160501_cutted_by_bounding_box.tif'
cut_by_bounding_box(531961.73, 535987.34, 4894164.57, 4888631.61)
print 'cutter.py script done!'
Best Answer
There are now Python modules easier to use for that, as rasterio
Therefore from Masking raster with a polygon feature in Rasterio Cookbook
Result
Alternatives for files
1) You can simply use the json or the geojson modules to import the geometry
2) you can use the ogr module with a shapefile
3) you can also use the Fiona module
New
You can use a filter as in the script of Luke in How to set a spatial filter with Python/GDAL?. Instead of cutting, you filter the input.
You only need to add the crs