Hi Even and community,
here's the link to the dataset that causes the error:
http://www.econ.brown.edu/students/Omer_Ozak_Munoz/Datasets/afgdist.zip Data
I was using the following code which used another dataset to create a mask,
but that's not the issue, so it could be simply deleted. I really appreciate
your help.
Best,
Omer
#############################
import os, sys, time
try:
from osgeo import gdal, gdalconst
from osgeo.gdalconst import *
import numpy
use_numeric = False
except:
import gdal, gdalconst
from gdalconst import *
import Numeric
use_numeric = True
# start timing
startTime=time.time()
pais = 'afg'
file= pais+'.txt'
file = open(file,'w')
# prepare gdal for importing ARCGIS binary grid
#gdal.AllRegister()
driver = gdal.GetDriverByName('AIG')
driver.Register()
if use_numeric:
# set the directory
dire='/direct/'
else:
dire='/Volumes/direct/'
# the file to open
fn = dire+'timehours'
mask = gdal.Open(fn, GA_ReadOnly)
if mask is None:
print('Could not open' + fn)
sys.exit(1)
else:
#determine info of grid
mcols = mask.RasterXSize
mrows = mask.RasterYSize
mbands = mask.RasterCount
mgeotransform = mask.GetGeoTransform()
moriginX = mgeotransform[0]
moriginY = mgeotransform[3]
mpixelWidth = mgeotransform[1]
mpixelHeight = mgeotransform[5]
mband = mask.GetRasterBand(1)
# set the directory
if use_numeric:
os.chdir('/direct')
else:
os.chdir('/Volumes/direct')
# the file to open
fn = pais+'dist'
ds = gdal.Open(fn, GA_ReadOnly)
if ds is None:
print('Could not open ' + fn)
sys.exit(1)
else:
#determine info of grid
cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount
geotransform = ds.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
band = ds.GetRasterBand(1)
# read the whole data
#data = band.ReadAsArray(0, 0, cols, rows)
# read only a couple of lines in (x,y)
x=-2684213
y=6146882
#x=1000000
#y=1000000
suma=0
count=0
maxi=0
xOffset = int((x-originX) / pixelWidth)
yOffset = int((y-originY) / pixelHeight)
y = -3883447
yend = int((y-originY) / pixelHeight)
xBsize = 100
yBsize = 100
#print 'mrows='+str(mrows)+' mcol='+str(mcols)
#print 'rows='+str(rows)+' col='+str(cols)
#print 'xoffset=' + str(xOffset) +', yoffset=' + str(yOffset)
for j in range(xOffset,cols,xBsize):
if j + xBsize <= cols:
ncols=xBsize
else:
ncols=cols-j
for i in range (yOffset,yend,yBsize):
if i+yBsize<=yend:
nrows=yBsize
else:
nrows=yend-i
if use_numeric:
data = band.ReadAsArray(j, i, ncols ,
nrows).astype(Numeric.Float)
mdata = mband.ReadAsArray(j, i, ncols , nrows)
mdata = Numeric.greater(data,0)
data = data * mdata
count = Numeric.sum(Numeric.sum(mdata))+count
suma = Numeric.sum(Numeric.sum(data))+suma
maxi = max(max(max(data)), maxi)
else:
data = band.ReadAsArray(j, i, ncols ,
nrows).astype(numpy.float32)
mdata = mband.ReadAsArray(j, i, ncols , nrows)
mdata = numpy.greater(mdata,0)
data = data * mdata
count = numpy.sum(numpy.sum(mdata))+count
suma = numpy.sum(numpy.sum(data))+suma
maxi=numpy.maximum(numpy.max(data),maxi)
#print 'i=' +str(i)+' j='+str(j)+' count=' +str(count)
print(' media=' + str(suma / count)+' max=' +str(maxi)+'
suma='+str(suma)+' count='+str(count))
file.write('code,mean,max,sum \n')
file.write(pais+ str(suma/count)+','+str(maxi)+','+str(suma)+' \n')
#print 'bands=' + str(bands) + ', originX=' + str(originX) + ',
originY=' + str(originY)
print('xoffset=' + str(xOffset) +', yoffset=' + str(yOffset))
endTime=time.time()
print('It took ' + str(startTime-endTime) + ' seconds')
file.close()
--
View this message in context:
http://n2.nabble.com/Re-corrupt-block-on-AIG-dataset-tp4568372p4571794.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
