473,289 Members | 1,917 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,289 software developers and data experts.

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.  
Sep 17 '21 #1
0 3543

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

Similar topics

5
by: DaveF | last post by:
IAny idea what this error is? C:\Documents and Settings\dfetrow\My Documents\Visual Studio Projects\FTPTimer\Service1.cs(15): The type or namespace name 'TimedItem' could not be found (are you...
3
by: trouling | last post by:
When trying to Create new ASP.NET project: ERROR: Visual Studio .NET has detected that the specified Web server is not running ASP.NET version 1.1 You will be unable to run ASP.NET web...
3
by: Regnab | last post by:
I'm not sure if this is some random bug in the system, but I keep getting the #Name? error on one of my forms when I try to link one field to another. For example, I have one field called...
1
by: 05l8kr | last post by:
I'm learning how to use Visual Studio 6.0 using beginners workshop, This program wants me to enter the following program and get it to compile and execute(run). Any help would be great. This is...
13
by: Protoman | last post by:
I'm getting an error: 10 C:\Dev-Cpp\Enigma.cpp no match for 'operator<' in 'i < (+cleartext)->std::basic_string<_CharT, _Traits, _Alloc>::end ()' Code: Enigma.hpp...
10
by: Durduran | last post by:
I can't seem to figure out the origin of this error and how I could go about fixing it -- the program that compiled on an older version of the OS no longer compiles with the upgraded version. ...
5
by: Alan Silver | last post by:
Hello, Server configuration: Windows 2003 Server SP2 SQL Server 2000 SP4 ..NET v2.0.50727 just built up a new server using the same configuration as my current one. I even used the same CDs...
1
by: CapCity | last post by:
I just grabbed a ASP.Net 2.0 (C# code behind) application from the repository. Builds OK. When I try to run it, in either debug or release mode, I get a dialog that says: "Unable to start...
4
by: munkee | last post by:
As topic says. Under normal database usage when I click quit I write to a log the current date/time. This is just for audit purposes to keep track of when people enter/leave the database. However,...
2
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 7 Feb 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:30 (7.30PM). In this month's session, the creator of the excellent VBE...
0
by: MeoLessi9 | last post by:
I have VirtualBox installed on Windows 11 and now I would like to install Kali on a virtual machine. However, on the official website, I see two options: "Installer images" and "Virtual machines"....
0
by: Aftab Ahmad | last post by:
Hello Experts! I have written a code in MS Access for a cmd called "WhatsApp Message" to open WhatsApp using that very code but the problem is that it gives a popup message everytime I clicked on...
0
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
1
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
0
by: jfyes | last post by:
As a hardware engineer, after seeing that CEIWEI recently released a new tool for Modbus RTU Over TCP/UDP filtering and monitoring, I actively went to its official website to take a look. It turned...
0
by: ArrayDB | last post by:
The error message I've encountered is; ERROR:root:Error generating model response: exception: access violation writing 0x0000000000005140, which seems to be indicative of an access violation...
1
by: PapaRatzi | last post by:
Hello, I am teaching myself MS Access forms design and Visual Basic. I've created a table to capture a list of Top 30 singles and forms to capture new entries. The final step is a form (unbound)...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.