preselecting rasters based on area of interest

michaelcollins1

I have almost 1900 orthos (200m x 200m, in ecw format) sent to me for a project. Does anyone have a model, python script  or method to preselect only those ortho images that intersect the study area? Loading all of them onto a map will take forever, and creating a raster dataset isn't working for some reason. I only have a basic license

 

thanks in advance

XanderBakker

I guess you could take some of the code in the script below, whichs runs through a server and creates a featureclass with polygons of the extents of the rasters found. You could intersect the extent with the polygon of your AOI and create a list of the relevant rasters, or add them directly to an open ArcMap session...

import arcpy
import os

def main():
    from collections import OrderedDict
    arcpy.env.overwriteOutput = True

    # output featureclass
    fc_out = r"C:\Path\To\Your\output.gdb\featureclass"

    workspace  = r"\\server\folder\\with\Rasters"
    rasters = []

    # 'extent', 'spatialReference',
    props_ras = ['name', 'extension','path',
                'catalogPath','bandCount', 'height', 'width', 'meanCellHeight',
                'meanCellWidth', 'pixelType', 'noDataValue', 'isInteger',
                'dataElementType', 'dataType']

    dct_prop_fld = OrderedDict({'name': 'RasName', 'extension': 'RasExt','path': 'RasPath',
                    'catalogPath': 'RasCatPath','bandCount': 'BandCount',
                    'height': 'Height', 'width': 'Width',
                    'meanCellHeight': 'CellH', 'meanCellWidth': 'CellW',
                    'pixelType': 'PixelType', 'noDataValue': 'NoDataValue',
                    'isInteger': 'IsInteger', 'dataElementType': 'DataElementType',
                    'dataType': 'DataType'})

    fld_srname = 'SR_name'
    fld_srwkid = 'SR_wkid'
    flds_sr = [fld_srname, fld_srwkid]

    flds_corr = ['RasName', 'RasCatPath', 'RasExt', 'BandCount', 'RasPath']

    dct_fldtype = {'RasName': ('TEXT', 50), 'RasExt': ('TEXT', 10),
                  'RasPath': ('TEXT', 255), 'RasCatPath': ('TEXT', 255),
                  'BandCount': ('LONG', None), 'Height': ('LONG', None), 'Width': ('LONG', None),
                    'CellH': ('DOUBLE', None), 'CellW': ('DOUBLE', None), 'PixelType': ('TEXT', 20),
                    'NoDataValue': ('DOUBLE', None), 'IsInteger': ('TEXT', 20),
                    'DataElementType': ('TEXT', 100), 'DataType': ('TEXT', 100),
                    fld_srname: ('TEXT', 255), fld_srwkid: ('LONG', None)}

    # create empty output fc
    sr = arcpy.SpatialReference(4686) # change this
    fc_ws, fc_name = os.path.split(fc_out)
    arcpy.CreateFeatureclass_management(fc_ws, fc_name, "POLYGON", spatial_reference=sr)

    # add fields
    for fldname, tpl_type in dct_fldtype.items():
        addField(fc_out, fldname, tpl_type)

    # start insert cursor
    flds = ['SHAPE@']
    flds.extend(dct_fldtype.keys())
    cnt = 0
    with arcpy.da.InsertCursor(fc_out, flds) as curs:

        # walk rasters
        walk = arcpy.da.Walk(workspace, topdown=True, datatype="RasterDataset")
        for dirpath, dirnames, filenames in walk:
            for filename in filenames:
                ras = os.path.join(dirpath, filename)

                try:
                    lst_row = []
                    desc = arcpy.Describe(ras)
                    if hasattr(desc, "bandcount"):
                        cnt += 1
                        print "Processing raster {0}: {1}".format(cnt, ras)
                        band_cnt = desc.bandCount
                        if band_cnt > 1:
                            band = 1
                            ras_ok = ras
                            desc_ok = arcpy.Describe(ras_ok)
                            ras = os.path.join(ras, "Band_{0}".format(band))
                            desc = arcpy.Describe(ras)
                        else:
                            desc_ok = desc

                        # describe ras sr and extent
                        sr = getAttribute(desc, 'spatialReference')
                        if not sr is None:
                            sr_name = sr.name
                            sr_wkid = sr.factoryCode
                        else:
                            sr_name = 'Unknown'
                            sr_wkid = -1

                        ext = getAttribute(desc, 'extent')
                        if not ext is None:
                            # create polygon
                            polygon = createPolygonFromExtent(ext)
                            lst_row.append(polygon)

                            for fldname in flds[1:]:
                                if fldname == fld_srname:
                                    lst_row.append(sr_name)
                                elif fldname == fld_srwkid:
                                    lst_row.append(sr_wkid)
                                else:
                                    prop = getKeyByValue(dct_prop_fld, fldname)
                                    if prop == None:
                                        print "prop == None: " + fldname
                                    else:
                                        if fldname in flds_corr:
                                            val = getAttribute(desc_ok, prop)
                                        else:
                                            val = getAttribute(desc, prop)
                                        if not val is None:
                                            lst_row.append(val)
                                        else:
                                            lst_row.append(None)

                            row = tuple(lst_row)
                            curs.insertRow(row)

                except Exception as e:
                    print e
                    continue

def getKeyByValue(dct, search_val):
    key_found = None
    for key, val in dct.items():
        if val == search_val:
            key_found = key
            break
    return key_found

def addField(fc_out, fldname, tpl_type):
    # TEXT | FLOAT | DOUBLE | SHORT | LONG | DATE | BLOB | RASTER | GUID
    if len(arcpy.ListFields(fc_out, wild_card=fldname)) == 0:
        if tpl_type[0] == 'TEXT':
            arcpy.AddField_management(fc_out, fldname, 'TEXT', field_length=tpl_type[1])
        else:
            arcpy.AddField_management(fc_out, fldname, tpl_type[0])

def getAttribute(desc, prop):
    if hasattr(desc, prop):
        return eval("desc.{0}".format(prop))
    else:
        return None

def createPolygonFromExtent(ext):
    corners = ['lowerLeft', 'upperLeft', 'upperRight', 'lowerRight', 'lowerLeft']
    arr = arcpy.Array()
    for corner in corners:
        arr.add(eval("ext.{0}".format(corner)))
    return arcpy.Polygon(arr, ext.spatialReference)

if __name__ == '__main__':
    main()

Apart from this, you could create an unmanaged Raster catalog I think...