需要帮助使用 gdal/osgeo 在 python 中计算区域统计数据

如何解决需要帮助使用 gdal/osgeo 在 python 中计算区域统计数据

嘿伙计们，我真的需要一些帮助来弄清楚为什么我在尝试计算 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]),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[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))
        else:
            zstats.append(setFeatureStats( \
                id,\
                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 插件，我不会收到任何错误。那么我做错了什么？我真的很感激这方面的一些帮助，这让我发疯。 谢谢