Hi everyone,
I have a list of rasters named: "ih1", "ih2", "ih3", ..., "ih28" that represent the interpolated rainfall height in the corresponding time step (hour). I want to calculate the cumulative rainfall height in each time step which is the sum of the rainfall in the same step and the cumulative rain of the previous one, e.g. CumRain(i) = ih(i) + CumRain(i-1). I am trying to code this in arcpy yet the "i" thing does not work for rasters.
If you have any idea how this could be solved, let me know.
Thanks a lot !
Convert the raster s to numpy arrays
RasterToNumPyArray—ArcGIS Pro | Documentation
a0
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
a1
array([[16, 17, 18, 19],
[20, 21, 22, 23],
[24, 25, 26, 27],
[28, 29, 30, 31]])
a2
array([[32, 33, 34, 35],
[36, 37, 38, 39],
[40, 41, 42, 43],
[44, 45, 46, 47]])
np.cumsum((a0, a1, a2), axis=0)
array([[[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]],
[[16, 18, 20, 22],
[24, 26, 28, 30],
[32, 34, 36, 38],
[40, 42, 44, 46]],
[[48, 51, 54, 57],
[60, 63, 66, 69],
[72, 75, 78, 81],
[84, 87, 90, 93]]], dtype=int32)
# ---- or np.nancumsum(a, axis=0) if you have nodata cells
Now if you want an intermediate time step, just slice it
csum = np.nancumsum(a, axis=0)
csum[1]
array([[16, 18, 20, 22],
[24, 26, 28, 30],
[32, 34, 36, 38],
[40, 42, 44, 46]], dtype=int32)
Then if you need raster s back
NumPyArrayToRaster—ArcGIS Pro | Documentation
Thank you, I try the code below:
from arcpy import env
from arcpy import env
from arcpy.sa import *
import os
import numpy
arcpy.env.workspace = "C:/MSc_thesis/modf"
idw = arcpy.ListRasters("ih*", "GRID")
for i in idw:
inras = Raster(i)
lowerLeft = arcpy.Point(inras.extent.XMin,inras.extent.YMin)
cellSize = inras.meanCellWidth
arr = arcpy.RasterToNumPyArray(inras,lowerLeft,cellSize, nodata_to_value=0)
csum = numpy.cumsum(arr)
newras = arcpy.NumPyArrayToRaster(csum,lowerLeft,cellSize,value_to_nodata=0)
but it returns the following error:
ValueError: Argument in_array: A two or three dimensional NumPy array is required.
you aren't accumulating the rasters, hence you get one value
arrays_in = []
for i in idw:
inras = Raster(i)
lowerLeft = arcpy.Point(inras.extent.XMin,inras.extent.YMin)
cellSize = inras.meanCellWidth
arr = arcpy.RasterToNumPyArray(inras,lowerLeft,cellSize, nodata_to_value=0)
arrays_in.append(arr)
main array = np.array(arrays_in)
c_sum = np.nancumsum(main_array, axis=0)
for i in c_sum:
newras = arcpy.NumPyArrayToRaster(*** you need to provide an output name for "i"
.You have to get all the arrays together, then if you need individual output rasters, then you need to provide a name for each output raster.
Also, don't forget the "axis" for the cumulative sum.
numpy.nancumsum returns error in my arcmap version so I replaced it with numpy.cumsum. Also, I might have an error in the code because now an ERROR 999998: Unexpected Error is appeared. I am not really familiar with all this stuff, sorry for the repetitive replies.
arrays_in = []
for i in idw:
inras = Raster(i)
lowerLeft = arcpy.Point(inras.extent.XMin,inras.extent.YMin)
cellSize = inras.meanCellWidth
arr = arcpy.RasterToNumPyArray(inras,lowerLeft,cellSize, nodata_to_value=0)
arrays_in.append(arr)
main_arr = numpy.array(arrays_in)
c_sum = numpy.cumsum(main_arr, axis = 0)
for j in c_sum:
newras = arcpy.NumPyArrayToRaster(j,lowerLeft,cellSize, value_to_nodata=0 )
newras.save(os.path.join(arcpy.env.workspace,"c" + str(j) + ".tif"))
what is one of the paths? sounds like it is bailing there. workspace would have to be a folder since you can't save tiffs in a gdb
Hi Cami,
If your hourly rasters are named sequentially as ih1.tif, ih2.tif, ih3.tif, etc. and you want to create cumulative rainfall rasters CumRain1.tif, CumRain2.tif, CumRain3.tif, etc. where conceptually:
CumRain1.tif = ih1.tif
CumRain2.tif = CumRain1.tif + ih2.tif
CumRain3.tif = CumRain2.tif + ih3.tif
and so on...
You can do it using a few simple steps:
1. Copy the first hourly rainfall raster to the first cumulative rainfall raster
2. loop through from 2nd hourly rainfall raster to the last hourly rainfall raster
3. for each hour create a cumulative rainfall raster by adding the hourly rainfall raster to the previous cumulative rainfall raster.
Please see the code snippet below.
Please note you can also use the Cell Statistics tool to add many rasters using the SUM option.
Hope it helps. Please let me know if you have any questions.
Thanks
Noman
#Code Snippet
import os
import sys
import arcpy
from arcpy.sa import *
#Input, Scratch, and Output Workspace
current_path = os.path.dirname(__file__)
data_folder = os.path.join(current_path, 'data')
out_folder = os.path.join(current_path, 'output')
scratch_folder = os.path.join(current_path, 'scratch')
#Geoprocessing environmnet
arcpy.env.scratchWorkspace = scratch_folder
arcpy.env.workspace = out_folder
arcpy.env.overwriteOutput = True
#Checkout Spatial Analyst extenson
arcpy.CheckOutExtension("spatial")
#Copy 1st hourly rainfall to 1st cumulative rainfall.
firstHourlyRainfallRasterName = "ih"+str(1)+".tif"
firstHourlyRainfallRaster = os.path.join(data_folder, firstHourlyRainfallRasterName)
firstCumulativeRainfallRasterName = "CumRain"+str(1)+".tif"
firstCumulativeRainfallRaster = os.path.join(out_folder, firstCumulativeRainfallRasterName)
arcpy.management.CopyRaster(firstHourlyRainfallRaster, firstCumulativeRainfallRaster)
#Loop through 2nd to last (10th) hourly rainfall.
for i in range(2, 11):
#Specifying input and output raster
PreviousCumulativeRainfallRasterName = "CumRain"+str(i-1)+".tif"
PreviousCumulativeRainfallRaster = os.path.join(out_folder, PreviousCumulativeRainfallRasterName)
hourlyRainfallRasterName = "ih"+str(i)+".tif"
hourlyRainfallRaster = os.path.join(data_folder, hourlyRainfallRasterName)
cumulativeRainfallRasterName = "CumRain"+str(i)+".tif"
cumulativeRainfallRaster = os.path.join(out_folder, cumulativeRainfallRasterName)
#adding hourly rainfall raster with previous cummulative rainfall raster
CumRain = Raster(PreviousCumulativeRainfallRaster) + Raster(hourlyRainfallRaster)
CumRain.save(cumulativeRainfallRaster)