[gdal-dev] Re: corrupt block on AIG dataset

[ozak]

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.