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

[Even Rouault]

Omer,

I've downloaded the dataset and run gdalinfo -checksum on it successfully (I'm 
using gdal-trunk - equivalent to 1.7 as far as AIG is concerned on 64bit 
Linux)

Now I'm running your script (with a slight modification to use the dataset as 
mask too), and until here (takes quite a long time to run!), I see no 
problem.

Not sure why you get errors...

Le Sunday 14 February 2010 22:41:40 ozak, vous avez écrit :
> 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()