469,282 Members | 1,650 Online
Bytes | Developer Community
New Post

Home Posts Topics Members FAQ

Post your question to a community of 469,282 developers. It's quick & easy.

I'm trying to figure out why I always get this error while trying to calculate zonal statistics in 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.
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?

Expand|Select|Wrap|Line Numbers
  1. import gdal
  2. from osgeo import ogr
  3. import os
  4. import numpy as np
  5. import csv
  6. import sys
  7.  
  8. if os.name == 'nt':
  9.     path=os.path.dirname(sys.argv[0])
  10. if os.name == 'posix':
  11.     path=os.path.dirname(os.path.abspath(__file__))
  12.  
  13. def boundingBoxToOffsets(bbox, geot):
  14.     col1 = int((bbox[0] - geot[0]) / geot[1])
  15.     col2 = int((bbox[1] - geot[0]) / geot[1]) + 1
  16.     row1 = int((bbox[3] - geot[3]) / geot[5])
  17.     row2 = int((bbox[2] - geot[3]) / geot[5]) + 1
  18.     return [row1, row2, col1, col2]
  19.  
  20.  
  21. def geotFromOffsets(row_offset, col_offset, geot):
  22.     new_geot = [
  23.         geot[0] + (col_offset * geot[1]),
  24.         geot[1],
  25.         0.0,
  26.         geot[3] + (row_offset * geot[5]),
  27.         0.0,
  28.         geot[5]
  29.     ]
  30.     return new_geot
  31.  
  32.  
  33. def setFeatureStats(fid, min, max, mean, median, sd, sum, count, names=["min", "max", "mean", "median", "sd", "sum", "count", "id"]):
  34.     featstats = {
  35.         names[0]: min,
  36.         names[1]: max,
  37.         names[2]: mean,
  38.         names[3]: median,
  39.         names[4]: sd,
  40.         names[5]: sum,
  41.         names[6]: count,
  42.         names[7]: fid,
  43.     }
  44.     return featstats
  45.  
  46. mem_driver = ogr.GetDriverByName("Memory")
  47. mem_driver_gdal = gdal.GetDriverByName("MEM")
  48. shp_name = "temp"
  49.  
  50. # fn_raster = "C:/cadastru2/task1B/T35TNK_20210512T085601_NDVI.tif"
  51. fn_raster =os.path.join(path,'ndvi','T35TNK_20210512T085601_NDVI.tif')
  52. fn_zones = "D:/cadastru2/task3/agri terenuri/Agri_Terenuri_NDVI.shp"
  53.  
  54.  
  55. r_ds = gdal.Open(fn_raster,1)
  56. p_ds = ogr.Open(fn_zones)
  57.  
  58. lyr = p_ds.GetLayer()
  59. geot = r_ds.GetGeoTransform()
  60. nodata = r_ds.GetRasterBand(1).GetNoDataValue()
  61.  
  62. zstats = []
  63.  
  64. p_feat = lyr.GetNextFeature()
  65. niter = 0
  66.  
  67. while p_feat:
  68.     if p_feat.GetGeometryRef() is not None:
  69.         if os.path.exists(shp_name):
  70.             mem_driver.DeleteDataSource(shp_name)
  71.         tp_ds = mem_driver.CreateDataSource(shp_name)
  72.         tp_lyr = tp_ds.CreateLayer('polygons', None, ogr.wkbPolygon)
  73.         tp_lyr.CreateFeature(p_feat.Clone())
  74.         offsets = boundingBoxToOffsets(p_feat.GetGeometryRef().GetEnvelope(), \
  75.                                        geot)
  76.         new_geot = geotFromOffsets(offsets[0], offsets[2], geot)
  77.  
  78.         tr_ds = mem_driver_gdal.Create( \
  79.             "", \
  80.             offsets[3] - offsets[2], \
  81.             offsets[1] - offsets[0], \
  82.             1, \
  83.             gdal.GDT_Byte)
  84.  
  85.         tr_ds.SetGeoTransform(new_geot)
  86.         gdal.RasterizeLayer(tr_ds, [1], tp_lyr, burn_values=[1])
  87.         tr_array = tr_ds.ReadAsArray()
  88.  
  89.         r_array = r_ds.GetRasterBand(1).ReadAsArray( \
  90.             offsets[2], \
  91.             offsets[0], \
  92.             offsets[3] - offsets[2], \
  93.             offsets[1] - offsets[0])
  94.  
  95.         id = p_feat.GetFID()
  96.  
  97.         if r_array is not None:
  98.             maskarray = np.ma.MaskedArray( \
  99.                 r_array, \
  100.                 maskarray=np.logical_or(r_array==nodata, np.logical_not(tr_array)))
  101.  
  102.             if maskarray is not None:
  103.                 zstats.append(setFeatureStats( \
  104.                     id, \
  105.                     maskarray.min(), \
  106.                     maskarray.max(), \
  107.                     maskarray.mean(), \
  108.                     np.ma.median(maskarray), \
  109.                     maskarray.std(), \
  110.                     maskarray.sum(), \
  111.                     maskarray.count()))
  112.             else:
  113.                 zstats.append(setFeatureStats( \
  114.                     id, \
  115.                     nodata, \
  116.                     nodata, \
  117.                     nodata, \
  118.                     nodata, \
  119.                     nodata, \
  120.                     nodata, \
  121.                     nodata))
  122.         else:
  123.             zstats.append(setFeatureStats( \
  124.                 id, \
  125.                 nodata, \
  126.                 nodata, \
  127.                 nodata, \
  128.                 nodata, \
  129.                 nodata, \
  130.                 nodata, \
  131.                 nodata))
  132.  
  133.         tp_ds = None
  134.         tp_lyr = None
  135.         tr_ds = None
  136.  
  137.         p_feat = lyr.GetNextFeature()
  138.  
  139. # fn_csv = "C:/cadastru2/task1B/zstats.csv"
  140. fn_csv = os.path.join(path,'zstats.csv')
  141. col_names = zstats[0].keys()
  142. with open(fn_csv, 'w', newline='') as csvfile:
  143.     writer = csv.DictWriter(csvfile, col_names)
  144.     writer.writeheader()
  145.     writer.writerows(zstats)
  146.  
  147.  
  148.  
  149.  
  150.  
4 Weeks Ago #1
0 2809

Post your reply

Sign in to post your reply or Sign up for a free account.

Similar topics

5 posts views Thread by DaveF | last post: by
3 posts views Thread by trouling | last post: by
13 posts views Thread by Protoman | last post: by
1 post views Thread by CapCity | last post: by
1 post views Thread by CARIGAR | last post: by
reply views Thread by zhoujie | last post: by
By using this site, you agree to our Privacy Policy and Terms of Use.