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