Exportación TIFF por lotes, QGIS

Question

Intento exportar 30 TIFF diferentes a partir de un TIFF y 30 capas de cuadros delimitadores. utilizando QGIS.

¿Quizás alguien tenga una pista para mí?

Edición1: Tengo un mapa en formato TIFF y necesito exportar 30 diferentes áreas pequeñas de la misma en formato TIFF. Me cansé de uno por uno desde el menú Archivo, que funciona, pero lo necesito en formato por lotes. Intenté buscar alguna información relevante pero no tuve éxito.

Answer 1

Esto requiere una herramienta de geoprocesamiento personalizada: Todas las características de un shapefile se recortan del raster de origen y se guardan como un archivo independiente con las mismas propiedades. Puede ejecutar este script en el editor qgis python después de cambiar la ruta y los nombres de archivo.

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