I have a shapefile with many polygons and many raster datasets (tif). Now, for each polygon I want to extract its area from the each raster datasets overlapping the polygon and save it to a new file (tif).
The code below does the job. However, since not many polygons and rasters overlay I create a lot of "empty" new tifs. Thus, I'm searching for a way to check whether a certain polygon overlays (a part) of a certain raster before running Extract By Mask.
Sidenote: I use arcpy.env.extent = "MAXOF" since otherwise the non-overlapping cases create an "ExecuteError: ERROR 010092: Invalid output extent" during Extract By Mask
import arcpy from arcpy import env import os if arcpy.CheckExtension("Spatial") == "Available": arcpy.CheckOutExtension("Spatial") else: # raise a custom exception raise LicenseError arcpy.CheckOutExtension("3D") arcpy.env.extent = "MAXOF" # path name variables output = "OUTPUTFOLDER" if not os.path.exists(output): os.makedirs(output) fc="POLYGON_INPUT_FC" crusor = arcpy.SearchCursor(fc) fieldName = "Plot_ID" env.workspace = "WORKSPACE_FOLDER" rasterList = arcpy.ListRasters("*", "tif") print "rasterList created" for poly in crusor: mask = poly.Plot_ID # Create feature layer of current mask polygon print mask whereClause = '"' + fieldName + '" = '+ str(mask) print whereClause arcpy.MakeFeatureLayer_management(fc, 'currentMask', whereClause) for raster in rasterList: rastername=raster[:-4] # Create and save the clipped raster outRaster = arcpy.sa.ExtractByMask(raster, 'currentMask') outRaster.save(output + "\\" + str(mask) +"_" + rastername + ".tif") arcpy.Delete_management('currentMask')
أكثر...
The code below does the job. However, since not many polygons and rasters overlay I create a lot of "empty" new tifs. Thus, I'm searching for a way to check whether a certain polygon overlays (a part) of a certain raster before running Extract By Mask.
Sidenote: I use arcpy.env.extent = "MAXOF" since otherwise the non-overlapping cases create an "ExecuteError: ERROR 010092: Invalid output extent" during Extract By Mask
import arcpy from arcpy import env import os if arcpy.CheckExtension("Spatial") == "Available": arcpy.CheckOutExtension("Spatial") else: # raise a custom exception raise LicenseError arcpy.CheckOutExtension("3D") arcpy.env.extent = "MAXOF" # path name variables output = "OUTPUTFOLDER" if not os.path.exists(output): os.makedirs(output) fc="POLYGON_INPUT_FC" crusor = arcpy.SearchCursor(fc) fieldName = "Plot_ID" env.workspace = "WORKSPACE_FOLDER" rasterList = arcpy.ListRasters("*", "tif") print "rasterList created" for poly in crusor: mask = poly.Plot_ID # Create feature layer of current mask polygon print mask whereClause = '"' + fieldName + '" = '+ str(mask) print whereClause arcpy.MakeFeatureLayer_management(fc, 'currentMask', whereClause) for raster in rasterList: rastername=raster[:-4] # Create and save the clipped raster outRaster = arcpy.sa.ExtractByMask(raster, 'currentMask') outRaster.save(output + "\\" + str(mask) +"_" + rastername + ".tif") arcpy.Delete_management('currentMask')
أكثر...