Index: wflow-py/Scripts/read_arcinfo_files.py =================================================================== diff -u -r9dd1a78b1a2e6aa0e67e970d6235e4923b8bfc42 -r325bb5b03289359b9ad2c8b814e66c0bd0d53dd7 --- wflow-py/Scripts/read_arcinfo_files.py (.../read_arcinfo_files.py) (revision 9dd1a78b1a2e6aa0e67e970d6235e4923b8bfc42) +++ wflow-py/Scripts/read_arcinfo_files.py (.../read_arcinfo_files.py) (revision 325bb5b03289359b9ad2c8b814e66c0bd0d53dd7) @@ -3,7 +3,7 @@ # Created on August 17, 2000 # by Keith Cherkauer # -# This python script contains library functions whic read and write +# This python script contains library functions whic read and write # standard Arc/Info ASCII grid files, returning a data dictionary # with information from the file. # @@ -15,11 +15,13 @@ from string import atoi from string import atof -def read_ARCINFO_ASCII_grid(filename, FLAG="UNFILTERED", INTflag=0, - Extent={ "North": 0, - "South": 0, - "East": 0, - "West": 0}): + +def read_ARCINFO_ASCII_grid( + filename, + FLAG="UNFILTERED", + INTflag=0, + Extent={"North": 0, "South": 0, "East": 0, "West": 0}, +): """This routine reads a standard Arc/Info ASCII grid file and returns the coordinates and values in a data dictionary. An optional flag passed to the routine determines if values for @@ -29,41 +31,50 @@ the boundaries of a smaller sub-grid that will be extracted.""" try: - f = open(filename,"r") + f = open(filename, "r") except IOError, E: - print "ERROR: Unable to open or read the ArcInfo grid file %s" % filename + print "ERROR: Unable to open or read the ArcInfo grid file %s" % filename print E - fileTable = { "Ncells" : 0 } - return ( fileTable ) + fileTable = {"Ncells": 0} + return fileTable fileTable = { - "ncols" : atoi(split(f.readline())[1]), - "nrows" : atoi(split(f.readline())[1]), - "xllcorner" : atof(split(f.readline())[1]), - "yllcorner" : atof(split(f.readline())[1]), - "cellsize" : atof(split(f.readline())[1]), - "NODATA_value" : atoi(split(f.readline())[1]), - } + "ncols": atoi(split(f.readline())[1]), + "nrows": atoi(split(f.readline())[1]), + "xllcorner": atof(split(f.readline())[1]), + "yllcorner": atof(split(f.readline())[1]), + "cellsize": atof(split(f.readline())[1]), + "NODATA_value": atoi(split(f.readline())[1]), + } print fileTable if Extent["North"] and Extent["South"] and Extent["East"] and Extent["West"]: # redefine grid size for desired extent - nrows = int(( Extent["North"] - Extent["South"] ) / fileTable["cellsize"]) - ncols = int(( Extent["East"] - Extent["West"] ) / fileTable["cellsize"]) + nrows = int((Extent["North"] - Extent["South"]) / fileTable["cellsize"]) + ncols = int((Extent["East"] - Extent["West"]) / fileTable["cellsize"]) else: # define full extent - Extent["North"] = fileTable["yllcorner"] + fileTable["cellsize"]*fileTable["nrows"] + Extent["North"] = ( + fileTable["yllcorner"] + fileTable["cellsize"] * fileTable["nrows"] + ) Extent["South"] = fileTable["yllcorner"] - Extent["East"] = fileTable["xllcorner"] + fileTable["cellsize"]*fileTable["ncols"] - Extent["West"] = fileTable["xllcorner"] + Extent["East"] = ( + fileTable["xllcorner"] + fileTable["cellsize"] * fileTable["ncols"] + ) + Extent["West"] = fileTable["xllcorner"] ncols = fileTable["ncols"] nrows = fileTable["nrows"] # compute number of cells required fileTable["Ncells"] = nrows * ncols # allocate memory for all cells - fileTable["cell"] = [ fileTable["NODATA_value"] ]*fileTable["Ncells"] + fileTable["cell"] = [fileTable["NODATA_value"]] * fileTable["Ncells"] - print "North %f, South %f, East %f, West %f" % ( fileTable["yllcorner"] + fileTable["cellsize"]*fileTable["nrows"], fileTable["yllcorner"], fileTable["xllcorner"] + fileTable["cellsize"]*fileTable["ncols"], fileTable["xllcorner"] ) + print "North %f, South %f, East %f, West %f" % ( + fileTable["yllcorner"] + fileTable["cellsize"] * fileTable["nrows"], + fileTable["yllcorner"], + fileTable["xllcorner"] + fileTable["cellsize"] * fileTable["ncols"], + fileTable["xllcorner"], + ) print Extent print fileTable["nrows"], nrows print fileTable["ncols"], ncols @@ -75,71 +86,88 @@ # get current line from file line = f.readline() # compute current latitude - lat = (fileTable["yllcorner"] + ( fileTable["nrows"] - i ) * fileTable["cellsize"] - fileTable["cellsize"] / 2.) + lat = ( + fileTable["yllcorner"] + + (fileTable["nrows"] - i) * fileTable["cellsize"] + - fileTable["cellsize"] / 2. + ) if lat >= Extent["South"] and lat <= Extent["North"]: tmprcnt = tmprcnt + 1 # if latitude falls within defined extent, split lines to get column values - if ( INTflag ): tmpvals = map(atoi,split(line)) - else: tmpvals = map(atof,split(line)) - if ( len(tmpvals) != fileTable["ncols"] ): - print "ERROR: Number of items in row %i (%i) does not match the number of defined columns (%i)." % (i, len(tmpvals), fileTable["ncols"]) - for j in range(fileTable["ncols"]): + if INTflag: + tmpvals = map(atoi, split(line)) + else: + tmpvals = map(atof, split(line)) + if len(tmpvals) != fileTable["ncols"]: + print "ERROR: Number of items in row %i (%i) does not match the number of defined columns (%i)." % ( + i, + len(tmpvals), + fileTable["ncols"], + ) + for j in range(fileTable["ncols"]): # compute longitude for current cell - lng = (fileTable["xllcorner"] + ( j ) * fileTable["cellsize"] + fileTable["cellsize"] / 2.) + lng = ( + fileTable["xllcorner"] + + (j) * fileTable["cellsize"] + + fileTable["cellsize"] / 2. + ) if lng >= Extent["West"] and lng <= Extent["East"]: - if tmprcnt == 1: tmpccnt = tmpccnt + 1 + if tmprcnt == 1: + tmpccnt = tmpccnt + 1 # if longitude within extent boundaries, store current location try: fileTable["cell"][cellidx] = { - "lat" : lat, - "lng" : lng, - "value" : tmpvals[j] - } + "lat": lat, + "lng": lng, + "value": tmpvals[j], + } cellidx = cellidx + 1 except IndexError, errstr: # did not allocate enough memory, add additional cells - fileTable["cell"] = fileTable["cell"] + [ { - "lat" : lat, - "lng" : lng, - "value" : tmpvals[j] - } ] + fileTable["cell"] = fileTable["cell"] + [ + {"lat": lat, "lng": lng, "value": tmpvals[j]} + ] cellidx = cellidx + 1 - del(line) + del (line) - print "Number of rows filled: %i of %i" % ( tmprcnt, nrows ) - print "Number of cols filled: %i of %i" % ( tmpccnt, ncols ) - print "Number of cells filled: %i of %i" % ( cellidx, fileTable["Ncells"] ) - if tmprcnt != nrows: nrows = tmprcnt - if tmpccnt != ncols: ncols = tmpccnt + print "Number of rows filled: %i of %i" % (tmprcnt, nrows) + print "Number of cols filled: %i of %i" % (tmpccnt, ncols) + print "Number of cells filled: %i of %i" % (cellidx, fileTable["Ncells"]) + if tmprcnt != nrows: + nrows = tmprcnt + if tmpccnt != ncols: + ncols = tmpccnt if cellidx < fileTable["Ncells"]: fileTable["cell"] = fileTable["cell"][:cellidx] - if cellidx != fileTable["Ncells"]: fileTable["Ncells"] = cellidx - + if cellidx != fileTable["Ncells"]: + fileTable["Ncells"] = cellidx + if FLAG == "FILTERED": - if fileTable["NODATA_value"] == 0: - fileTable["cell"] = filter(lambda x: x["value"] != 0, fileTable["cell"]) - else: - fileTable["cell"] = filter(lambda x: x["value"] != -9999, fileTable["cell"]) - fileTable["Ncells"] = len(fileTable["cell"]) + if fileTable["NODATA_value"] == 0: + fileTable["cell"] = filter(lambda x: x["value"] != 0, fileTable["cell"]) + else: + fileTable["cell"] = filter(lambda x: x["value"] != -9999, fileTable["cell"]) + fileTable["Ncells"] = len(fileTable["cell"]) - # reset grid boundaries to agree with defined extent + # reset grid boundaries to agree with defined extent fileTable["ncols"] = ncols fileTable["nrows"] = nrows fileTable["xllcorner"] = Extent["West"] fileTable["yllcorner"] = Extent["South"] f.close() - return ( fileTable ) + return fileTable + def write_ARCINFO_ASCII_grid(filename, gridTable, INTflag=0): """This routine writes a standard Arc/Info ASCII grid file values of which are stored in a data dictionary.""" try: - f = open(filename,"w") + f = open(filename, "w") except IOError, E: - print "ERROR: Unable to open or write the ArcInfo grid file %s" % filename - return ( 0 ) + print "ERROR: Unable to open or write the ArcInfo grid file %s" % filename + return 0 f.write("ncols\t%i\n" % gridTable["ncols"]) f.write("nrows\t%i\n" % gridTable["nrows"]) @@ -151,24 +179,43 @@ idx = 0 CellsWritten = 0 for i in range(gridTable["nrows"]): - lat = (gridTable["yllcorner"] + ( gridTable["nrows"] - i ) * gridTable["cellsize"] - gridTable["cellsize"] / 2.) + lat = ( + gridTable["yllcorner"] + + (gridTable["nrows"] - i) * gridTable["cellsize"] + - gridTable["cellsize"] / 2. + ) tmpstr = "" - for j in range(gridTable["ncols"]): - lng = (gridTable["xllcorner"] + ( j ) * gridTable["cellsize"] + gridTable["cellsize"] / 2.) - if idx < gridTable["Ncells"] and (abs(lat-gridTable["cell"][idx]["lat"])<=gridTable["cellsize"]/2. and abs(lng-gridTable["cell"][idx]["lng"])<=gridTable["cellsize"]/2.): - if ( INTflag or gridTable["cell"][idx]["value"] == gridTable["NODATA_value"]): + for j in range(gridTable["ncols"]): + lng = ( + gridTable["xllcorner"] + + (j) * gridTable["cellsize"] + + gridTable["cellsize"] / 2. + ) + if idx < gridTable["Ncells"] and ( + abs(lat - gridTable["cell"][idx]["lat"]) <= gridTable["cellsize"] / 2. + and abs(lng - gridTable["cell"][idx]["lng"]) + <= gridTable["cellsize"] / 2. + ): + if ( + INTflag + or gridTable["cell"][idx]["value"] == gridTable["NODATA_value"] + ): tmpstr = tmpstr + "%i " % gridTable["cell"][idx]["value"] - else: tmpstr = tmpstr + "%f " % gridTable["cell"][idx]["value"] + else: + tmpstr = tmpstr + "%f " % gridTable["cell"][idx]["value"] if gridTable["cell"][idx]["value"] != gridTable["NODATA_value"]: CellsWritten = CellsWritten + 1 idx = idx + 1 else: tmpstr = tmpstr + "%i " % gridTable["NODATA_value"] - f.write(tmpstr[:-1]+"\n") + f.write(tmpstr[:-1] + "\n") f.close() - print "%s: %i cells of %i contain data." % ( filename, CellsWritten, - gridTable["ncols"]*gridTable["nrows"] ) + print "%s: %i cells of %i contain data." % ( + filename, + CellsWritten, + gridTable["ncols"] * gridTable["nrows"], + ) - return (1) + return 1