I have 30 rasters with 0 and 1 values and I want to combine them into a single raster using the combine tool. However, combine tool has a 20 raster limit, so my here is my current approach and code:

  1. Split the rasters into two batches of 15 rasters and run the combine
  2. Next, run combine on the outputs created in step #1 to create final output.

        import arcpy,os,time
        from arcpy import env
        from arcpy.sa import *
        # Main Directories
        CCDC = "D:/CCDC/"
        CCDC_bi = "D:/LC/CCDC_bi/"
        years = ['1985', '1986', '1987', '1988', '1989', '1990', 
            '1991', '1992', '1993', '1994', '1995', '1996', 
            '1997', '1998', '1999', '2000', '2001', '2002', 
            '2003', '2004', '2005', '2006', '2007', '2008', 
            '2009', '2010', '2011', '2012', '2013', '2014', 
        for year in years:
            output = os.path.join(CCDC_bi, "year_" + year + "_fix.tif")
            output1 = os.path.join(CCDC_bi, "year_" + year + "_bi.tif")
            if arcpy.Exists(output) and arcpy.Exists(output1):
                print output + " // " + output1 + " >>> Already Exist! "
                print "Binary Raster Creation Started..."
                outCon4 = Con(IsNull(output), 0, 1)
                arcpy.CopyRaster_management(outCon4, output1,"","","","NONE","NONE","4_BIT","NONE","NONE","TIFF")
                print "Binary Raster Creation Complete"
        comb_list_1 = [] #Placeholder for all combine paths
        comb_list_2 = [] #Placeholder for all combine paths
        bi_years_1 = years[3:-14]
        bi_years_2 = years[-14:-1]
        for i in bi_years_1:
            inRas = os.path.join(CCDC_bi, "year_" + i + "_bi.tif")
            if arcpy.Exists(inRas):
                comb_list_1.append(CCDC_bi + "year_" + i + "_bi.tif")
        for i in bi_years_2:
            inRas = os.path.join(CCDC_bi, "year_" + i + "_bi.tif")
            if arcpy.Exists(inRas):
                comb_list_2.append(CCDC_bi + "year_" + i + "_bi.tif")
        comb1 = ";".join(comb_list_1)
        comb2 = ";".join(comb_list_2)
        start_time = time.time()
        print 'Starting Combine - 1'
        outCombine_1 = Combine(comb1)
        print("--- Combine 1 Complete Complete %s seconds ---" % (time.time() - start_time))
        start_time = time.time()
        print 'Starting Combine - 2'
        outCombine_2 = Combine(comb2)
        print("--- Combine 2 Complete Complete %s seconds ---" % (time.time() - start_time))
        start_time = time.time()
        print 'Starting Combine Final'
        finalCombine = Combine("D:/LC/CCDC_bi/combine1.tif;D:/LC/CCDC_bi/combine2.tif")
        print("--- Combine Final Complete Complete %s seconds ---" % (time.time() - start_time))

I am not sure if this approach will work, I'm currently running this model now to look at the final output. Alternatively, I have tried running combine on these rasters on R but it runs out of memory. Has anyone tried combining large numbers of rasters and if so what method or strategy worked best for them?

  • I think you should let your code run and only worry about seeking an alternative if it fails to complete. You will then be in a position to tell us not just what you want to do and what you have tried, but also where you are actually stuck. – PolyGeo Dec 14 '16 at 22:05
  • How are you going to use output please? I am asking because even combination of 5 values is hard to analyse, I cannot imagine dealing with 30. Number of unique combinations is enormous – FelixIP Dec 14 '16 at 22:18
  • @FelixIP, the rasters represent binary annual land use change. The plan is to create an output that has land change values for all years in a single raster where we can use the columns for past-future prediction analysis. It allows us to keep our data in a tabular and spatial format throughout our analysis as we explore different scenarios, so when we have to create change rasters from the output we don't have to worry about joining tables to rasters. In fact, your suggestion on my recent question about looping through ~16 mil. records was a product of 13 rasters combined! – cptpython Dec 15 '16 at 1:11

Something between this lines (untested !):

for year in range(1985,2016):
  output = arcpy.Raster(os.path.join(CCDC_bi, "year_" + str(year) + "_fix.tif"))
  binary = Con(IsNull(output), 0, Power(2,year-1985))
  toSave = os.path.join(CCDC_bi, "year_" + str(year) + "_fixSaved.tif"))
  listOfRasters.append(toSave) in
outCellStats = CellStatistics(listOfRasters, "SUM", "DATA")

will do. When you'll get down to analysis, you'll need to think "binary' e.g. value of 9 in result would mean that changes happened in 1985 and 1988, because 9 converted to binary is 1001.

|improve this answer|||||
  • you'll need to fix your 3rd line. It will fail because you're concatenating a string and an int. You need to add str(year). Same with the fifth line. – Fezter Dec 15 '16 at 1:43
  • If I understand this correctly, you convert all continous rasters to binary rasters then calculate sum of all the cells using cell statistics. The end product should be a single column of binary numbers which indicate the year of first change. This might work but I would need to check with the statistical modeler, because he created binary rasters i.e. "output1" using focal statistics and region group. – cptpython Dec 15 '16 at 3:14
  • No. I convert your 0,1 rasters to something else. First year raster will contain 0,1 second 0,2 fourth 0,8 tenth 0,1024 – FelixIP Dec 15 '16 at 3:28
  • Ah, I see. I will report on this soon. – cptpython Dec 15 '16 at 3:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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