If you have the 3D Analyst extension, use the Surface Volume tool.
It's possible to automate this calculation in Python to get volumes at different heights (from an incremental height change). (Note: this example is with version 9.3, but it isn't too difficult to convert it to arcpy version 10).
# ArcGIS version 9.3
import arcgisscripting
import re
# Edit this section
step_size = 10 # metres; reduce this if you want finer resolution
num_steps = 20 # increase this, if you need more steps
workspace = r'C:\Some\path\to\folder'
raster_file = 'mytopo.tif'
gp = arcgisscripting.create()
gp.CheckOutExtension ('3D')
gp.workspace = workspace
# Start from the minimum elevation
min_rast = gp.GetRasterProperties(raster_file, 'MINIMUM')
# Start CSV output of elevation/volume relationship
print('step, height, elev, volume')
for step in range(num_steps):
height = float(step*step_size)
elev = min_rast + height
gp.SurfaceVolume_3d(raster_file, '', 'BELOW', elev)
result = gp.GetMessages()
volume = float(re.findall(r'Volume= *([\d\.]+)', result)[0])
print(', '.join([repr(x) for x in [step + 1, height, elev, volume]]))
produces this CSV output:
step, height, elev, volume
1, 0.0, 195.72373962402344, 0.0
2, 10.0, 205.72373962402344, 89383.136320076999
3, 20.0, 215.72373962402344, 537710.29799682996
4, 30.0, 225.72373962402344, 1616009.8535998999
5, 40.0, 235.72373962402344, 4022313.9784607999
6, 50.0, 245.72373962402344, 10139496.423829
7, 60.0, 255.72373962402344, 22423139.656699002
8, 70.0, 265.72373962402344, 125592027.17357001
9, 80.0, 275.72373962402344, 601472924.93704998
10, 90.0, 285.72373962402344, 1454118895.066
11, 100.0, 295.72373962402344, 2659321774.0138998
12, 110.0, 305.72373962402344, 4271651250.0858998
13, 120.0, 315.72373962402344, 6145982306.6198997
14, 130.0, 325.72373962402344, 8248103003.7566004
15, 140.0, 335.72373962402344, 10540914755.195
16, 150.0, 345.72373962402344, 13094286080.875
17, 160.0, 355.72373962402344, 15855580342.691999
18, 170.0, 365.72373962402344, 18894132090.799
19, 180.0, 375.72373962402344, 22168466487.139
20, 190.0, 385.72373962402344, 25702467387.048
Following @Dango's idea I created and tested (on small rasters with the same extent and cell size) the following code:
import arcpy, numpy
inRaster = r"C:\tmp\RastersArray.gdb\InRaster"
inRaster2 = r"C:\tmp\RastersArray.gdb\InRaster2"
##Get properties of the input raster
inRasterDesc = arcpy.Describe(inRaster)
#coordinates of the lower left corner
rasXmin = inRasterDesc.Extent.Xmin
rasYmin = inRasterDesc.Extent.Ymin
# Cell size, raster size
rasMeanCellHeight = inRasterDesc.MeanCellHeight
rasMeanCellWidth = inRasterDesc.MeanCellWidth
rasHeight = inRasterDesc.Height
rasWidth = inRasterDesc.Width
##Calculate coordinates basing on raster properties
#create numpy array of coordinates of cell centroids
def rasCentrX(rasHeight, rasWidth):
coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)
return coordX
inRasterCoordX = numpy.fromfunction(rasCentrX, (rasHeight,rasWidth)) #numpy array of X coord
def rasCentrY(rasHeight, rasWidth):
coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)
return coordY
inRasterCoordY = numpy.fromfunction(rasCentrY, (rasHeight,rasWidth)) #numpy array of Y coord
#combine arrays of coordinates (although array for Y is before X, dstack produces [X, Y] pairs)
inRasterCoordinates = numpy.dstack((inRasterCoordY,inRasterCoordX))
##Raster conversion to NumPy Array
#create NumPy array from input rasters
inRasterArrayTopLeft = arcpy.RasterToNumPyArray(inRaster)
inRasterArrayTopLeft2 = arcpy.RasterToNumPyArray(inRaster2)
#flip array upside down - then lower left corner cells has the same index as cells in coordinates array
inRasterArray = numpy.flipud(inRasterArrayTopLeft)
inRasterArray2 = numpy.flipud(inRasterArrayTopLeft2)
# combine coordinates and value
inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))
#add values from second raster
rasterValuesArray = numpy.dstack((inRasterFullArray, inRasterArray2.T))
Based on @hmfly code, you can have access to desired values:
(height, width, dim )=rasterValuesArray.shape
for row in range(0,height):
for col in range(0,width):
#now you have access to single array of values for one cell location
Unfortunately there's one 'but' - the code is right for NumPy arrays which can be handled by system memory. For my system (8GB), the largest array was about 9000,9000.
As my experience doesn't let me provide more help, you can consider some suggestions about dealing wiht large arrays:
https://stackoverflow.com/questions/1053928/python-numpy-very-large-matrices
arcpy.RasterToNumPyArray method allows to specify the subset of raster converted to NumPy array (ArcGIS10 help page) what can be useful when chunking large dataset into submatrices.
Best Answer
I don't know of a specific tool to perform the conversion you are looking to do, but I would guess that it should be relatively straightforward to calculate (please check my math).
Q in mm = 0.001 Q in m. Now you can multiply Q in m by the area to get your volume:0.001m * 10,000m2 = 10m3per secondand you needper day, your conversion factor is60*60*24 = 86,400 seconds per day.(10m3/sec)*86,400 sec/day = 864,000 m3/dayThis can all be done using the raster calculator tool.