Select to view content in your preferred language

Implementing merge() in ArcPy map algebra

5368
1
09-01-2014 09:17 AM
curtvprice
MVP Esteemed Contributor

The merge function (which we've used for years in GRID, ArcView, and ArcGIS 8x, 9x) simply doesn't exist in ArcPy map algebra, unless I'm missing something...

Here's my hack. Note this is set up so it should work inside a map algebra expression. It also returns an integer raster if all the inputs are integer, otherwise float.

Just asking the community (and especially Esri SA team folks - I know you're out there) whether this is a viable approach.

import os

import arcpy

def MergeRasters(raster_list, out_raster=None): 

    """ArcPy Map Algebra implementation of the old merge() function.

     Curtis Price, cprice@usgs.gov

    """ 

   

    # output raster name

    if not out_raster:

        out_raster = arcpy.CreateScratchName(

            "mrg", "", "raster", arcpy.env.workspace) 

       

    # determine raster type for output

    ptype = "32_BIT_SIGNED" 

    for ras in raster_list: 

        if not arcpy.Describe(ras).isInteger: 

            ptype = "32_BIT_FLOAT" 

            break 

       

    # create temporary raster

    tmpRaster = arcpy.CreateScratchName("", "mrg", "raster", 

                                          arcpy.env.scratchFolder) 

    arcpy.CreateRasterDataset_management(os.path.dirname(tmpRaster), 

                              os.path.basename(tmpRaster), 

                              pixel_type=ptype) 

    arcpy.Mosaic_management(raster_list, tmpRaster, "FIRST") 

   

    # Copy the raster, enforcing environment (extent, snap, mask, etc) 

    mRaster = ApplyEnvironment(tmpRaster)

    mRaster.save(out_raster)

   

    # Clean up 

    arcpy.Delete_management(tmpRaster) 

    return mRaster

Note, this may work sort of like the merge() function, but is not true arcpy map algebra local operator, so we lose out on the benefits of using it with complex statements like

Int(Merge(raster1, SetNull(raster2, 2), raster3 + 2))

Please vote up my idea: Bring back Spatial Analyst Merge tool

UPDATES

curtvprice - script bug fix

curtvprice - Used ApplyEnvironment function  instead of just multiplying by one

curtvprice - added note about complex statements and link to ideas.esri.com

Tags (2)
0 Kudos
1 Reply
curtvprice
MVP Esteemed Contributor

Here's a new implementation. I think this is a much better replacement for merge().

(I did get some help on this from one of the Spatial Analyst developer luminaries.)

import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")

def Merge(raster_list): 
    """Merge a list of rasters. The first non-NoData value found (cell by cell) is kept.
       Merge([raster1, raster2, raster3, ...])

    This is another implementation of ArcGIS 9x Map Algebra merge() function"""
    if type(raster_list) == str:
        raster_list = raster_list.split(";")
    if len(raster_list) < 2:
        raise Exception("More than one raster required")
    rlist = [Raster(arcpy.Describe(ras).catalogPath) for ras in raster_list]
    merge_ras = rlist[0]
    for ras in rlist[1:]:
        merge_ras = Con(IsNull(merge_ras), ras, merge_ras)
    return merge_ras