1

I have raster format radar data which contains the depth of rainfall in each grid cell. A radar image is taken every 5 minutes, but I want to know the total depth of rainfall in each day.

The radar data filenames are identical apart from the date and time. So for example:

metoffice-c-band-rain-radar_uk_201805312140_5km-composite.gif

For ease, I have separated the radar data by month into different folders (so one folder for May, one folder for June etc.)

I'm using Python to solve this problem. My code is designed to take the characters from the filename that indicate day, and store it as [day], I then create a list with all filenames for the same day being grouped together. I then use the list for each day as an input to cellstatistics. However, Python crashes everytime I run this tool or it responds with an unexpected error 999998.

import os, sys, arcpy
from collections import defaultdict
from arcpy import env
from arcpy.sa import *
arcpy.env.cartographicCoordinateSystem = arcpy.SpatialReference("British National Grid")
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference("British National Grid")
arcpy.env.extent = r"C:\Users\bwj202\OneDrive - University of Exeter\GIS\Leakage\RuralTrunkMains.shp"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
env.workspace = os.getcwd()

dictdata = defaultdict(list)
for filename in os.listdir(os.getcwd()):
     if filename.endswith('.gif'):
          day = filename[37:39]
          dictdata[day].append(filename)      

Day1 = dictdata["01"]
Day2 = dictdata["02"]
Day3 = dictdata["03"]
Day4 = dictdata["04"]
Day5 = dictdata["05"]
Day6 = dictdata["06"]
Day7 = dictdata["07"]
Day8 = dictdata["08"]
Day9 = dictdata["09"]
Day10 = dictdata["10"]
Day11 = dictdata["11"]
Day12 = dictdata["12"]
Day13 = dictdata["13"]
Day14 = dictdata["14"]
Day15 = dictdata["15"]
Day16 = dictdata["16"]
Day17 = dictdata["17"]
Day18 = dictdata["18"]
Day19 = dictdata["19"]
Day20 = dictdata["20"]
Day21 = dictdata["21"]
Day22 = dictdata["22"]
Day23 = dictdata["23"]
Day24 = dictdata["24"]
Day25 = dictdata["25"]
Day26 = dictdata["26"]
Day27 = dictdata["27"]
Day28 = dictdata["28"]
Day29 = dictdata["29"]
Day30 = dictdata["30"]
Day31 = dictdata["31"]

#print(Day1)
#print(Day2)
#print(dictdata)

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day1], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day1.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day2], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day2.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day3], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day3.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day4], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day4.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day5], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day5.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day6], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day6.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day7], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day7.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day8], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day8.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day9], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day9.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day10], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day10.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day11], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day11.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day12], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day12.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day13], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day13.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day14], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day14.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day15], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day15.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day16], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day16.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day17], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day17.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day18], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day18.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day19], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day19.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day20], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day20.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day21], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day21.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day22], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day22.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day23], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day23.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day24], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day24.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day25], "SUM", "NODATA")

# Save the output 
outCellStatistics1.save("Day25.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day26], "SUM", "NODATA")

outCellStatistics1.save("Day26.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day27], "SUM", "NODATA")

outCellStatistics1.save("Day27.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day28], "SUM", "NODATA")

outCellStatistics1.save("Day28.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day29], "SUM", "NODATA")

outCellStatistics1.save("Day29.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day30], "SUM", "NODATA")

outCellStatistics1.save("Day30.tif")

# Execute CellStatistics
outCellStatistics1 = CellStatistics([Day31], "SUM", "NODATA")

outCellStatistics1.save("Day31.tif")

#exit()

I've also tried running this as a loop, but this generates an error that states that arcpy has tried to overwrite the same file (e.g. it writes Day2.tif and then tries to write Day2.tif a second time). Here is the alternative code:

import os, sys, arcpy
from collections import defaultdict
from arcpy import env
from arcpy.sa import *
arcpy.env.cartographicCoordinateSystem = arcpy.SpatialReference("British National Grid")
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference("British National Grid")
arcpy.env.extent = r"C:\Users\bwj202\OneDrive - University of Exeter\GIS\Leakage\RuralTrunkMains.shp"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
env.workspace = os.getcwd()

dictdata = defaultdict(list)
for filename in os.listdir(os.getcwd()):
     if filename.endswith('.asc'):
          day = filename[37:39]
          dictdata[day].append(filename)      
for k in dictdata:
  save = os.path.join(os.getcwd(), k + ".tif")
  outCellStatistics = CellStatistics([dictdata[k]], "SUM", "NODATA")
  outCellStatistics.save(save)
Improve this question
1
  • Don't know if this is a copy n paste error but the 3 lines after for k in dictdata: need to be indented.
    – Hornbydd
    Commented Nov 22, 2018 at 17:06

2 Answers 2

Reset to default
0

I think it may be the way you create the dictionary and populate it:

Replace this:

dictdata = defaultdict(list)

with

dictdata = {}

Replace this line 

dictdata[day].append(filename)

with

dictdata[day] = filename

Also the way your code is structured you feed exactly 1 raster into the CellStatistics tool and the sum of itself is itself... If you are trying to sum rasters from day 1 through to day 31 you need to be adding them to a LIST not a dictionary because then you feed in a LIST of rasters to the CellStatistics tool which is what the help file is stating under the syntax section.

Improve this answer
6
  • I see. Well with the code you suggest I get an error: Day10 = dictdata["10"] KeyError:"10".
    – Ben_123
    Commented Nov 22, 2018 at 17:32
  • If I instead use dictdata = [] (rather than {} so I create a list rather than a dictionary) I get the error: list indices must be integers, not str.
    – Ben_123
    Commented Nov 22, 2018 at 17:38
  • Incidentally, if I print(Day1), I receive a list of filenames that do indeed relate to the first Day. But when I use Day1 in a cellstatistics statement, I receive: RuntimeError: ERROR 999998: Unexpected Error.
    – Ben_123
    Commented Nov 22, 2018 at 18:02
  • I think you need to add a screen shot of your files to your question. Not sure why you would get a key error unless there are no day 10?
    – Hornbydd
    Commented Nov 22, 2018 at 18:54
  • I'd love to provide you with a screenshot...but there are 8,300 files in the directory. But it does indeed appear that there is no day 10, so I removed reference to day 10 in the code. I also inserted: print(Day1) and it does indeed receive only one reference.
    – Ben_123
    Commented Nov 23, 2018 at 10:48
0

So I can answer my own question here, largely with help from colleagues. The main problem was that there were memory issues (which was odd as I was able to run cellstatistics in ArcMap) but I also needed the full path. Some of the other changes are implemented purely for the code to be "cleaner". Here is the code:

import os, sys, arcpy
from collections import defaultdict
from arcpy import env
from arcpy.sa import *
arcpy.env.cartographicCoordinateSystem = arcpy.SpatialReference("British National Grid")
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference("British National Grid")
arcpy.env.extent = r"C:\Users\bwj202\OneDrive - University of Exeter\GIS\Leakage\RuralTrunkMains.shp"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
Data_folder = os.getcwd() + "/Data"
env.workspace = os.getcwd()
print (Data_folder)

dictdata = defaultdict(list)
for filename in os.listdir(Data_folder):
     if filename.endswith('.gif'):
          day = filename[37:39]
          dictdata[day].append(Data_folder + "/" + filename)

# print (dictdata)

Export_folder = os.getcwd() + "/Exports"

if os.path.exists(Export_folder):
    pass
else:
    arcpy.CreateFolder_management(os.getcwd(), "Exports")


print ("Start loop")

for key, value in dictdata.iteritems():
    label = str(key)
    # print value
    outCellStatistics = CellStatistics(value, "SUM", "DATA") 

    # # Save the output
    outCellStatistics.save(Export_folder + "/Day" + label + ".tif")


print ("finished")
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.