需要使用gdal/osgeo在python中计算分区统计信息的帮助吗

Question

嘿，伙计们，我真的需要一些帮助来弄清楚为什么我在试图计算gdal中的分区统计数据时总是出现这个错误：

ERROR 5: D:/cadastru2/task1B\ndvi\T35TNK_20210512T085601_NDVI.tif, band 1: Access window out of range in RasterIO().  Requested (20979,445042) of size 45x13 on raster of 10980x10980.

首先，我的代码是：

import gdal
from osgeo import ogr
import os
import numpy as np
import csv
import sys

if os.name == 'nt':
    path=os.path.dirname(sys.argv[0])
if os.name == 'posix':
    path=os.path.dirname(os.path.abspath(__file__))

def boundingBoxToOffsets(bbox, geot):
    col1 = int((bbox[0] - geot[0]) / geot[1])
    col2 = int((bbox[1] - geot[0]) / geot[1]) + 1
    row1 = int((bbox[3] - geot[3]) / geot[5])
    row2 = int((bbox[2] - geot[3]) / geot[5]) + 1
    return [row1, row2, col1, col2]


def geotFromOffsets(row_offset, col_offset, geot):
    new_geot = [
        geot[0] + (col_offset * geot[1]),
        geot[1],
        0.0,
        geot[3] + (row_offset * geot[5]),
        0.0,
        geot[5]
    ]
    return new_geot


def setFeatureStats(fid, min, max, mean, median, sd, sum, count, names=["min", "max", "mean", "median", "sd", "sum", "count", "id"]):
    featstats = {
        names[0]: min,
        names[1]: max,
        names[2]: mean,
        names[3]: median,
        names[4]: sd,
        names[5]: sum,
        names[6]: count,
        names[7]: fid,
    }
    return featstats

mem_driver = ogr.GetDriverByName("Memory")
mem_driver_gdal = gdal.GetDriverByName("MEM")
shp_name = "temp"

# fn_raster = "C:/cadastru2/task1B/T35TNK_20210512T085601_NDVI.tif"
fn_raster =os.path.join(path,'ndvi','T35TNK_20210512T085601_NDVI.tif')
fn_zones = "D:/cadastru2/task3/agri terenuri/Agri_Terenuri_NDVI.shp"


r_ds = gdal.Open(fn_raster,1)
p_ds = ogr.Open(fn_zones)

lyr = p_ds.GetLayer()
geot = r_ds.GetGeoTransform()
nodata = r_ds.GetRasterBand(1).GetNoDataValue()

zstats = []

p_feat = lyr.GetNextFeature()
niter = 0

while p_feat:
    if p_feat.GetGeometryRef() is not None:
        if os.path.exists(shp_name):
            mem_driver.DeleteDataSource(shp_name)
        tp_ds = mem_driver.CreateDataSource(shp_name)
        tp_lyr = tp_ds.CreateLayer('polygons', None, ogr.wkbPolygon)
        tp_lyr.CreateFeature(p_feat.Clone())
        offsets = boundingBoxToOffsets(p_feat.GetGeometryRef().GetEnvelope(), \
                                       geot)
        new_geot = geotFromOffsets(offsets[0], offsets[2], geot)

        tr_ds = mem_driver_gdal.Create( \
            "", \
            offsets[3] - offsets[2], \
            offsets[1] - offsets[0], \
            1, \
            gdal.GDT_Byte)

        tr_ds.SetGeoTransform(new_geot)
        gdal.RasterizeLayer(tr_ds, [1], tp_lyr, burn_values=[1])
        tr_array = tr_ds.ReadAsArray()

        r_array = r_ds.GetRasterBand(1).ReadAsArray( \
            offsets[2], \
            offsets[0], \
            offsets[3] - offsets[2], \
            offsets[1] - offsets[0])

        id = p_feat.GetFID()

        if r_array is not None:
            maskarray = np.ma.MaskedArray( \
                r_array, \
                maskarray=np.logical_or(r_array==nodata, np.logical_not(tr_array)))

            if maskarray is not None:
                zstats.append(setFeatureStats( \
                    id, \
                    maskarray.min(), \
                    maskarray.max(), \
                    maskarray.mean(), \
                    np.ma.median(maskarray), \
                    maskarray.std(), \
                    maskarray.sum(), \
                    maskarray.count()))
            else:
                zstats.append(setFeatureStats( \
                    id, \
                    nodata, \
                    nodata, \
                    nodata, \
                    nodata, \
                    nodata, \
                    nodata, \
                    nodata))
        else:
            zstats.append(setFeatureStats( \
                id, \
                nodata, \
                nodata, \
                nodata, \
                nodata, \
                nodata, \
                nodata, \
                nodata))

        tp_ds = None
        tp_lyr = None
        tr_ds = None

        p_feat = lyr.GetNextFeature()

# fn_csv = "C:/cadastru2/task1B/zstats.csv"
fn_csv = os.path.join(path,'zstats.csv')
col_names = zstats[0].keys()
with open(fn_csv, 'w', newline='') as csvfile:
    writer = csv.DictWriter(csvfile, col_names)
    writer.writeheader()
    writer.writerows(zstats)

问题是，如果我将光栅添加到Qgis中，并使用zonalStats插件，我不会得到任何错误。那么我做错了什么？我真的很感激你能帮我，这让我发疯了。 多谢各位