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
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...
Thanks, I'll take a look at the code, it's probably the way to go.
This is as far as I got with a model. I'm not sure what to add at this point, whereby if Top < x and Bottom > y .etc. then move raster to folder Z, or simply collect the name of the raster if it qualifies and save the list as a text file.