Sep-17-2021, 09:00 AM
I'm trying to figure out why I always get this error while trying to calculate zonal
Error: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.
The problem is if I add my raster to QGIS and use the zonalStats plugin I get no errors whatsoever. So what am I doing wrong?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)