Python Forum
I'm trying to figure out why I always get this error while trying to calculate zonal
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
I'm trying to figure out why I always get this error while trying to calculate zonal
#1
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)
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  I cannot figure out this error ErnestTBass 8 3,048 Aug-20-2020, 08:06 PM
Last Post: ErnestTBass
  I cannot figure out this error ErnestTBass 8 3,424 Apr-24-2020, 06:03 PM
Last Post: ErnestTBass
  I cannot figure out this error ErnestTBass 6 2,544 Mar-28-2020, 04:58 AM
Last Post: SheeppOSU
  Python error. Can't figure it out. ignatius 9 4,726 Oct-21-2018, 11:47 PM
Last Post: ignatius
  Can't figure out the syntax error, relatively simple code maho686868 3 3,083 Jul-08-2018, 03:43 PM
Last Post: volcano63

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020