5

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.

Original image Cutted by bounding box

Cutted by GeoJSON

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!'
  • I can not use the function: bbox, lat_lon_to_utm and explode. What libraries do I have to do this? – user84253 Oct 12 '16 at 9:27
  • You can use geojson module. First install it pip install geojson then import and try it import geojson as gj and gj.utils.coords() You' ll get list of coordinates like: [[38.5711773250159, 45.6353245353988], [38.5653971566629, 45.6384032362039], ... [38.5711773250159, 45.6353245353988]] For more info see an examples – Sant Oct 12 '16 at 9:31
7

There are now Python modules easier to use for that, as rasterio

Rasterio employs GDAL to read and writes files using GeoTIFF and many other formats. Its API uses familiar Python and SciPy interfaces and idioms like context managers, iterators, and ndarrays.

Therefore from Masking raster with a polygon feature in Rasterio Cookbook

enter image description here

import rasterio
from rasterio.mask import mask
# the polygon GeoJSON geometry
geoms = [{'type': 'Polygon', 'coordinates': [[(250204.0, 141868.0), (250942.0, 141868.0), (250942.0, 141208.0), (250204.0, 141208.0), (250204.0, 141868.0)]]}]
# load the raster, mask it by the polygon and crop it
with rasterio.open("test.tif") as src:
    out_image, out_transform = mask(src, geoms, crop=True)
out_meta = src.meta.copy()

# save the resulting raster  
out_meta.update({"driver": "GTiff",
    "height": out_image.shape[1],
    "width": out_image.shape[2],
"transform": out_transform})

with rasterio.open("masked.tif", "w", **out_meta) as dest:
    dest.write(out_image)

Result

enter image description here

Alternatives for files

1) You can simply use the json or the geojson modules to import the geometry

with open('poly.json') as data_file:    
    geoms= json.load(data_file)
print geoms
{u'type': u'Polygon', u'coordinates': [[[250204.0, 141868.0], [250942.0, 141868.0], [250942.0, 141208.0], [250204.0, 141208.0], [250204.0, 141868.0]]]}

2) you can use the ogr module with a shapefile

from osgeo import ogr
reader = ogr.Open('poly.json')
layer = reader.GetLayer()
feature = layer.GetFeature(0)
geoms =json.loads(feature.ExportToJson())['geometry']
print geoms
{u'type': u'Polygon', u'coordinates': [[[250204.0, 141868.0], [250942.0, 141868.0], [250942.0, 141208.0], [250204.0, 141208.0], [250204.0, 141868.0]]]}

3) you can also use the Fiona module

It is a Python wrapper for vector data access functions from the OGR library

import fiona
with fiona.open("poly.shp") as shapefile:
    geoms = [feature["geometry"] for feature in shapefile]
print geoms
[{'type': 'Polygon', 'coordinates': [[(250204.0, 141868.0), (250942.0, 141868.0), (250942.0, 141208.0), (250204.0, 141208.0), (250204.0, 141868.0)]]}]

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.

from osgeo import gdal
xmin,ymin,xmax,ymax = [250204.0, 141208.0, 250942.0, 141868.0]
def map_to_pixel(mx,my,gt):
   #Convert from map to pixel coordinates.
   #Only works for geotransforms with no rotation.
   px = int((mx - gt[0]) / gt[1]) #x pixel
   py = int((my - gt[3]) / gt[5]) #y pixel
   return px, py

def extent_to_offset(xmin,ymin,xmax,ymax,gt):
  pxmin,pymin = map_to_pixel(xmin,ymin,gt)
  pxmax,pymax = map_to_pixel(xmax,ymax,gt)
  return pxmin,pymin,abs(pxmax-pxmin),abs(pymax-pymin)

ds=gdal.Open('test.tif')
gt = ds.GetGeoTransform()

bands = ds.RasterCount
band_list = []
#
# Read in bands and store all the data in bandList
#
for i in range(bands):
   band = ds.GetRasterBand(i+1) # 1-based index
   # apply filter to only read the data in the bounding box (xmin, ...)
   data = band.ReadAsArray(xoff, yoff, xsize, ysize)
   band_list.append(data)

 driver = gdal.GetDriverByName("GTiff")
 dst_dst = driver.Create('result.tif', xsize, ysize, 4, gdal.GDT_Byte)
 for j in range(bands):
    data = band_list[j]
    dst_dst.GetRasterBand(j+1).WriteArray(data)

 ....
 dst_dst = None

You only need to add the crs

  • Nice solutions, but I forgot to said that I can only use GDAL and numpy modules. I'll edit my question and title. Sry. – Sant Jun 24 '16 at 10:27
  • look at Clip a GeoTiff with Shapefile in the Python GDAL/OGR Cookbook – gene Jun 24 '16 at 10:38
  • and Clip a Raster using a Shapefile – gene Jun 24 '16 at 10:46
  • yep I think it's only one way to cut tiff files by geojson without using any console utilities and other modules. I was stuck for 4 days with this issue. And I don't found any other solutions. Thx – Sant Jun 24 '16 at 12:12
  • look above my solution in New – gene Jun 24 '16 at 13:45
3

Here is my own solution. It works for any number of bands, any types of geometry(e.g. multipolygon) and works with images any zones!

import geojson as gj
from osgeo import ogr, osr, gdal

# Enable GDAL/OGR exceptions
gdal.UseExceptions()


# GDAL & OGR memory drivers
GDAL_MEMORY_DRIVER = gdal.GetDriverByName('MEM')
OGR_MEMORY_DRIVER = ogr.GetDriverByName('Memory')


def cut_by_geojson(input_file, output_file, shape_geojson):

    # Get coords for bounding box
    x, y = zip(*gj.utils.coords(gj.loads(shape_geojson)))
    min_x, max_x, min_y, max_y = min(x), max(x), min(y), max(y)

    # Open original data as read only
    dataset = gdal.Open(input_file, gdal.GA_ReadOnly)

    bands = dataset.RasterCount

    # Getting georeference info
    transform = dataset.GetGeoTransform()
    projection = dataset.GetProjection()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = -transform[5]

    # Getting spatial reference of input raster
    srs = osr.SpatialReference()
    srs.ImportFromWkt(projection)

    # WGS84 projection reference
    OSR_WGS84_REF = osr.SpatialReference()
    OSR_WGS84_REF.ImportFromEPSG(4326)

    # OSR transformation
    wgs84_to_image_trasformation = osr.CoordinateTransformation(OSR_WGS84_REF,
                                                                srs)
    XYmin = wgs84_to_image_trasformation.TransformPoint(min_x, max_y)
    XYmax = wgs84_to_image_trasformation.TransformPoint(max_x, min_y)

    # Computing Point1(i1,j1), Point2(i2,j2)
    i1 = int((XYmin[0] - xOrigin) / pixelWidth)
    j1 = int((yOrigin - XYmin[1]) / pixelHeight)
    i2 = int((XYmax[0] - xOrigin) / pixelWidth)
    j2 = int((yOrigin - XYmax[1]) / pixelHeight)
    new_cols = i2 - i1 + 1
    new_rows = j2 - j1 + 1

    # New upper-left X,Y values
    new_x = xOrigin + i1 * pixelWidth
    new_y = yOrigin - j1 * pixelHeight
    new_transform = (new_x, transform[1], transform[2], new_y, transform[4],
                     transform[5])

    wkt_geom = ogr.CreateGeometryFromJson(str(shape_geojson))
    wkt_geom.Transform(wgs84_to_image_trasformation)

    target_ds = GDAL_MEMORY_DRIVER.Create('', new_cols, new_rows, 1,
                                          gdal.GDT_Byte)
    target_ds.SetGeoTransform(new_transform)
    target_ds.SetProjection(projection)

    # Create a memory layer to rasterize from.
    ogr_dataset = OGR_MEMORY_DRIVER.CreateDataSource('shapemask')
    ogr_layer = ogr_dataset.CreateLayer('shapemask', srs=srs)
    ogr_feature = ogr.Feature(ogr_layer.GetLayerDefn())
    ogr_feature.SetGeometryDirectly(ogr.Geometry(wkt=wkt_geom.ExportToWkt()))
    ogr_layer.CreateFeature(ogr_feature)

    gdal.RasterizeLayer(target_ds, [1], ogr_layer, burn_values=[1],
                        options=["ALL_TOUCHED=TRUE"])

    # Create output file
    driver = gdal.GetDriverByName('GTiff')
    outds = driver.Create(output_file, new_cols, new_rows, bands,
                          gdal.GDT_Float32)

    # Read in bands and store all the data in bandList
    mask_array = target_ds.GetRasterBand(1).ReadAsArray()
    band_list = []

    for i in range(bands):
        band_list.append(dataset.GetRasterBand(i + 1).ReadAsArray(i1, j1,
                         new_cols, new_rows))

    for j in range(bands):
        data = np.where(mask_array == 1, band_list[j], mask_array)
        outds.GetRasterBand(j + 1).SetNoDataValue(0)
        outds.GetRasterBand(j + 1).WriteArray(data)

    outds.SetProjection(projection)
    outds.SetGeoTransform(new_transform)

    target_ds = None
    dataset = None
    outds = None
    ogr_dataset = None
  • Where does field_shape_geojson come from in wkt_geom = ogr.CreateGeometryFromJson(str(field_shape_geojson))? Do you mean to use shape_geojson which is passed as input parameter? – Kamrul Dec 4 '18 at 6:52
  • 1
    @Kamrul sure, fixed) – Sant Dec 4 '18 at 21:33

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.