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() _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
