2

I try to export 30 different TIFFs from one TIFF and 30 layer of boundaries boxes using QGIS.

Maybe someone has a hint for me?

Edit1: I have a map in TIFF format and I need to export 30 different small areas from it in TIFF format. I tired one by one from File menu, that works but I need it in batch format. I tried to search some relevant information but I had no success.

  • Please be more specific. Explain, what you tried so far, what went wrong, what you expect to achieve exactly, etc. – Erik Nov 7 at 13:17
  • I tried single exports from file menu, it's working but I need to export using batch. I expect to receive 30 tiff files with my own boundary sets that lay in layers/layer/table – Sasha Nov 7 at 13:44
  • 1
    Any what is wrong with clip raster by layer mask? – Erik Nov 7 at 14:21
  • Are your boundaries represented as polygons in a single vector file? – Walshe_d Nov 9 at 19:30
  • Yes they're all in the same layer – Sasha Nov 9 at 20:13
1

This calls for a custom geoprocessing tool: Alle the features in a shapefile are clipped form the source raster and saved as a separate file with same properties. You can run this script in the qgis python editor after changing the path and filenames.

import os, sys
from osgeo import gdal, gdalnumeric, ogr, osr
import numpy as np

# config, all your data is in OUT_FOLDER
OUT_FOLDER= r"c:\tmp"
IN_SHAPE = r"tiles.shp"
IN_RAS =  r"raster.tif"

#open shapefile
source_shape = os.path.join(out_path, IN_SHAPE)
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(source_shape, 0)
layer = dataSource.GetLayer()

# open raster
raster_path = os.path.join(OUT_FOLDER, IN_RAS)
print raster_path

srcImage = gdal.Open(raster_path)

# create the spatial reference of source raster
srs = osr.SpatialReference()
srs.ImportFromWkt(srcImage.GetProjectionRef())

# geotransform info
gt = srcImage.GetGeoTransform()
rb = srcImage.GetRasterBand(1)
gdal_datatype = rb.DataType
src_ds_nd = rb.GetNoDataValue()

# Get raster georeference info
xOrigin, yOrigin, pixelWidth, pixelHeight = gt[0], gt[3], gt[1], gt[5]
cellsize = pixelWidth

for feature in layer:
    n = 0
    geom = feature.GetGeometryRef()
    name = feature.GetField("name")
    env = geom.GetEnvelope()
    min_x, min_y, max_x, max_y = env[0], env[2], env[1], env[3]

    minX = round(min_x / cellsize, 0) * cellsize
    maxX = round(max_x / cellsize, 0) * cellsize 
    minY = round(min_y / cellsize, 0) * cellsize
    maxY = round(max_y / cellsize) * cellsize 

    ncol = int((maxX - minX) / cellsize)
    nrow= int((maxY - minY) / cellsize)

    tmp_ds = gdal.GetDriverByName('MEM').Create('', ncol, nrow, 1, gdal.GDT_Byte)

    tmp_ds.SetGeoTransform((
        minX, cellsize, 0,
        maxY, 0, -cellsize,
    ))

    #  Create for target raster the same projection as for the value raster
    dsmem = driver.CreateDataSource("/vsimem/temporarely.shp")
    l_out_mem = dsmem.CreateLayer("memlayer", srs, ogr.wkbPolygon)
       # Add an ID field
    idField = ogr.FieldDefn("id", ogr.OFTInteger)
    l_out_mem.CreateField(idField)

    # Create the feature and set values
    featureDefn = l_out_mem.GetLayerDefn()
    out_feature = ogr.Feature(featureDefn)
    out_feature.SetField("id", 1)
    out_feature.SetGeometry(geom.Clone())
    l_out_mem.CreateFeature(out_feature)

 # burn value 1 to new raster    
    tmp_ds.SetProjection(srs.ExportToWkt())
    gdal.RasterizeLayer(tmp_ds, [1], l_out_mem, burn_values=[1])

    dsmem = None
    # clip same extent from srcImage
    cellX = int((minX - xOrigin) / pixelWidth)
    cellY = int((yOrigin - maxY) / -pixelHeight)

    nparr = rb.ReadAsArray(cellX, cellY, ncol, nrow)


     # mask nieuwe array met rasterized layer
    bandmask = tmp_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0, 0, ncol, nrow).astype(np.float)
    datamask[datamask == 0] = 'nan'

    # Mask zone of raster
    if (nparr is not None):
        zonedata = nparr * datamask
        zonedata[zonedata == src_ds_nd] = 'nan'
        n = np.count_nonzero(~np.isnan(zonedata))

    if (nparr is None or n < 1) :
        QMessageBox.information(None, "No data is found in feature: "+name,"")
        continue

        # negeer NoData
    zonedata = zonedata.astype(np.float)
    zonedata[zonedata == src_ds_nd] = 'nan'

        # write the clipped raster to fiel    
    target_ds = gdal.GetDriverByName('GTiff').Create(os.path.join(OUT_FOLDER, name+'.tif'), nkol, nrij, 1, gdal_datatype)
    target_ds.SetGeoTransform((minX, cellsize, 0, maxY, 0, -cellsize, ))

    # Create for target raster the same projection as for the value raster
    target_ds.SetProjection(srs.ExportToWkt())
    outBand = target_ds.GetRasterBand(1)

    # write the data
    outBand.WriteArray(zonedata, 0, 0)
    outBand.FlushCache()
    outBand.SetNoDataValue(-9999)
    target_ds = None
    outBand = None 

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.