问题描述
嘿伙计们,我真的需要一些帮助来弄清楚为什么我在尝试计算 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 插件,我不会收到任何错误。那么我做错了什么?我真的很感激这方面的一些帮助,这让我发疯。 谢谢
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)