1

我正在使用多边形（shapefile，*.shp）裁剪光栅图像（Erdas Imagine Images，*.img），下面的示例Clip a Raster using a Shapefile。按照示例我有错误：

Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
IndexError: too many indices

帖子评论里面也报这个错误，但是没有报解决办法。使用此链接数据可以下载示例中使用的数据。当我使用 16 位图像时，图像是 8 位的。错误的起源可能是这种差异（8bit vs 16 bit）

import osgeo.gdal
from osgeo import gdal, gdalnumeric, ogr, osr
import numpy as np

def world2Pixel(geoMatrix, x, y):
    """
    Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
    the pixel location of a geospatial coordinate
    (source http://www2.geog.ucl.ac.uk/~plewis/geogg122/vectorMask.html)
    geoMatrix
    [0] = top left x (x Origin)
    [1] = w-e pixel resolution (pixel Width)
    [2] = rotation, 0 if image is "north up"
    [3] = top left y (y Origin)
    [4] = rotation, 0 if image is "north up"
    [5] = n-s pixel resolution (pixel Height)

    """
    ulX = geoMatrix[0]
    ulY = geoMatrix[3]
    xDist = geoMatrix[1]
    yDist = geoMatrix[5]
    rtnX = geoMatrix[2]
    rtnY = geoMatrix[4]
    pixel = np.round((x - ulX) / xDist).astype(np.int)
    line = np.round((ulY - y) / xDist).astype(np.int)
    return (pixel, line)

inFile = 'C:\\myimage.img'
poly = "C:\\ClipArea.shp"

# Open the image as a read only image
ds = osgeo.gdal.Open(inFile,gdal.GA_ReadOnly)
# Check the ds (=dataset) has been successfully open
# otherwise exit the script with an error message.
if ds is None:
    raise SystemExit("The raster could not openned")
# Get image georeferencing information.
geoMatrix = ds.GetGeoTransform()
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
# Get image georeferencing information.
geoMatrix = ds.GetGeoTransform()
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
# open shapefile (= border of are of interest)
shp = osgeo.ogr.Open(poly)
if not len(shp) is 1:
    print "The shapefile with more than 1 record"
    sys.exit(-1)
source_shp = ogr.GetDriverByName("Memory").CopyDataSource(shp, "")
# Create an OGR layer from a boundary shapefile
source_layer = source_shp.GetLayer(0)
# Convert the layer extent to image pixel coordinates
minX, maxX, minY, maxY = source_layer.GetExtent()
ulX, ulY = world2Pixel(geoMatrix, minX, maxY)
lrX, lrY = world2Pixel(geoMatrix, maxX, minY)
# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)

# Load the source data as a gdalnumeric array
srcArray = gdalnumeric.LoadFile(inFile)
clip = srcArray[:, ulY:lrY, ulX:lrX]

**Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
IndexError: too many indices**
4

0 回答 0