在读取面文件位置栅格数据时，出现Access window out of range in RasterIO()的问题

Question

@如雾如电 ，@谛凌，@Mr.Winter`
我想读取shp位置的遥感影像的各波段数值时，出现Access window out of range in RasterIO()，我想知道应该在哪里修改才能解决

def RandomForestClassification(ClassifyRaster, TrainRaster, TrainShp, outfile, blockSize=0, treenum=100, max_depth=10):
    rds = gdal.Open(ClassifyRaster)
    transform = (rds.GetGeoTransform())
    lX = transform[0]  # 左上角点
    lY = transform[3]
    rX = transform[1]  # 分辨率
    rY = transform[5]
    width = rds.RasterXSize
    height = rds.RasterYSize
    bX = lX + rX * width  # 右下角点
    bY = lY + rY * height
    BandsCount = rds.RasterCount
    clf = tbrf.createClassifier(TrainRaster, TrainShp)

def createClassifier(inraster, inshp, field: str = "Id", treenum: int = 100):
    rasterspatial = gdal.Open(inraster)
    spatial2 = osr.SpatialReference()
    spatial2.ImportFromWkt(rasterspatial.GetProjectionRef())
    shpspatial = ogr.Open(inshp)
    layer = shpspatial.GetLayer(0)
    spatial1 = layer.GetSpatialRef()
 
    ct = osr.CreateCoordinateTransformation(spatial1, spatial2)
    oFeature = layer.GetNextFeature()
    # 下面开始遍历图层中的要素
    geom = oFeature.GetGeometryRef()
    if geom.GetGeometryType() == ogr.wkbPoint:
        return createClassifierByPoint(inraster, inshp)
    k = geom.GetGeometryType()
    if geom.GetGeometryType() != ogr.wkbPolygon:
        print("样本必须为单部件多边形")
        return False
    trainX = list()
    trainY = list()
    print("读取样本")
    while oFeature is not None:
        geom = oFeature.GetGeometryRef()
        # print(geom.GetGeometryName())
        # print(geom.GetGeometryCount())
        # print(geom.GetPointCount())
        wkt = geom.ExportToWkt()
        points = WKTToPoints(wkt)
        polygonPoints = []
        value = oFeature.GetField(field)
        for point in points:
            pC = ct.TransformPoint(point.X, point.Y, 0)
            polygonPoints.append(Point(pC[0], pC[1]))
            arr, width, height, BandsCount, leftX, upY = GetSubRaster(inraster, polygonPoints)
 
**函数运行到这里，转到另一个函数**
 
def GetSubRaster(inraster, polygonPoints: list):
    polygonPoints.append(polygonPoints[0])  # 面多边形坐标封闭
    print("当前多边形节点数量：" + str(len(polygonPoints)))
    minX = 10000000000000
    maxX = -minX
    minY = 100000000000000000
    maxY = -minY
 
    for point in polygonPoints:
        if point.X < minX: minX = point.X
        if point.X > maxX: maxX = point.X
        if point.Y < minY: minY = point.Y
        if point.Y > maxY: maxY = point.Y
    leftX = minX
    upY = maxY
    rightX = maxX
    bottomY = minY
 
    rds = gdal.Open(inraster)
 
    transform = (rds.GetGeoTransform())
    lX = transform[0]  # 左上角点
    lY = transform[3]
    rX = transform[1]  # 分辨率
    rY = transform[5]
 
    wpos = int((leftX - lX) / rX)
    hpos = int((upY - lY) / rY)
 
    width = int((rightX - leftX) / rX)
    height = int((bottomY - upY) / rY)
    BandsCount = rds.RasterCount
    arr = rds.ReadAsArray(wpos, hpos, width, height)
    fixX = list()
 
    nodatavalue = rds.GetRasterBand(1).GetNoDataValue()
    for i in range(height):
        if height > 200:
            print(f"多边形裁剪进度：{round(((i + 1) / height) * 100, 4)}%")
 
        Y = upY + i * rY + .00001
        pointsindex = list()
        for k in range(len(polygonPoints) - 1):
            point1 = polygonPoints[k]
            point2 = polygonPoints[k + 1]
            if (point1.Y >= Y and point2.Y <= Y) or (point1.Y <= Y and point2.Y >= Y):
                pointsindex.append(k)
        for j in range(width):
            count = 0
            for m in (pointsindex):
                point1 = polygonPoints[m]
                point2 = polygonPoints[m + 1]
                X = leftX + j * rX + .00001
                if point1.X == point2.X:
                    intersectX = point1.X
                    if intersectX > X: count += 1
                else:
                    k = (point2.Y - point1.Y) / (point2.X - point1.X)
                    if k == 0:
                        if X < point1.X or X < point2.X:
                            count += 1
                    else:
                        intersectX = (Y - point1.Y) / k + point1.X
                        if intersectX > X: count += 1
 
            if count % 2 == 0:
                if BandsCount > 1:
                    for bc in range(BandsCount):
                    arr[bc][i][j] = (nodatavalue)
                else:
                    arr[i][j] = -1
 
    return arr, width, height, BandsCount, leftX, upY

出现的问题
ERROR 5: 1.tif: Access window out of range in RasterIO(). Requested (-6552644,-6552482) of size 65x76 on raster of 7434x6520.
TypeError: 'NoneType' object is not subscriptable

Answer 1

size 65x76 这个shp的范围可能超出了影像范围， 你需要检查下你的shp是否靠影像的边缘太近，如果太近了建议选择偏内部一点的地方来取范围，假如你一定要取那个位置的数据推荐你先扩充图像的边界，扩充方式可以是只填充一定宽度的0，也可以是在边缘一定范围内镜像，具体参考这里https://blog.csdn.net/qq_20373723/article/details/123142344