0

I have a folder of rasters (.tif) and I want to determine which rasters intersect a polygon. The rasters are single band with a value of "1" while the rest of the extent has a no data value (coloured green in the images) of "-3.39999995214e+38". I am only interested in the raster if cells containing "1" are within my polygon. I am not interested in a raster if my polygon intersects the raster's no data area or it's extent. This picture indicates that my polygon intersects the extent of this raster; however, there are no cell with a value of "1" within my polygon only no data values. In this example, I am not interested in this raster.

enter image description here

In the second picture, I am interested in this raster as cell values of "1" are within my polygon.

enter image description here

Any suggestions?

I am using ArcGIS 10.3.1 and Python 2.7. Here is my code:

# Python Script to identify if a raster has data within a prepresentative AOI 
# polygon.

# Import system modules
import arcpy
import os, sys
#import ogr, gdal

# Set Geoprocessing environments
arcpy.env.parallelProcessingFactor = "85%"
arcpy.env.overwriteOutput = "True"

#--------------------------------------------------------------------------------------------------------------
def myUtil(msg):
    # This little utility just lets us print the same messages no matter how it is run.
    #arcpy.AddMessage(msg) # for use in ArcCatalog
    print msg # for use in Eclipse
    return()
#--------------------------------------------------------------------------------------------------------------

# Script arguments
inFolder = r"G:\Junk\Refugia" 

# Identify the Landuse Framework Region.
LufArea = "NSR"

# Select the appropriate AOI polygon.
AOI     = r"U:\Projects\LUF\RDRP\InComingData\Refugia\Refugia_data.gdb\NRS_LUF_Region_Buffer_EcoRefugia" # For EcoRegion Rasters

# Create the appropriate folder to store the selected rasters.
oFolder = inFolder + os.sep + LufArea + "_Selected_EcoRegions" # the folder to copy selected EcoRegions

# Set the Workspace environment
arcpy.env.workspace = inFolder
myUtil("\nWorkspace set to : " + inFolder)

# Create the storage folder for rasters that have data within our stated AOI.
myUtil("\nCreating the raster selection storage folder.")
#NewFldr = LufArea + "_Selected_Species" # For Species Rasters
NewFldr = LufArea + "_Selected_Ecoregions" # For EcoRegion Rasters
arcpy.CreateFolder_management(inFolder, NewFldr)

# Define Raster Projection
myUtil("\nDefining\Verifying Raster Projections and Checking for Data within the AOI ...\n")
#for ThisRas in arcpy.ListRasters("*", "TIF"): # For Species Rasters
for ThisRas in arcpy.ListRasters("*eco_curr*", "TIF"): # For EcoRegion Rasters

    if arcpy.Describe(ThisRas).spatialreference.name == "Ecoregions - Lambert_Conformal_Conic_2SP": # For EcoRegion Rasters
        myUtil("  Raster %s Projection is good." % (ThisRas))
    else:
        #arcpy.DefineProjection_management(in_dataset=ThisRas, coor_system="PROJCS['Refugia - Lambert Azimuthal Equal Area',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Lambert_Azimuthal_Equal_Area'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-100.0],PARAMETER['Latitude_Of_Origin',45.0],UNIT['Meter',1.0]]")
        arcpy.DefineProjection_management(in_dataset=ThisRas, coor_system="PROJCS['Ecoregions - Lambert_Conformal_Conic_2SP',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Lambert_Conformal_Conic'],PARAMETER['false_easting',0.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',-95.0],PARAMETER['standard_parallel_1',49.0],PARAMETER['standard_parallel_2',77.0],PARAMETER['latitude_of_origin',0.0],UNIT['Meter',1.0]]")
        myUtil("  Raster %s Projection has been defined." % (ThisRas))

    # Process: Calculate Statistics
    myUtil ("    Calculating Raster %s statistics for the AOI portion only." % (ThisRas))
    arcpy.CalculateStatistics_management(ThisRas, "1", "1", "", "OVERWRITE", AOI)

    # Process: Get Maximum Value from Raster Properties
    myUtil ("      Determining the maximum value for Raster %s to see if there is any data in the AOI portion." % (ThisRas))
    MaxValResult = arcpy.GetRasterProperties_management(ThisRas, "MAXIMUM", "")
    MaxVal = MaxValResult.getOutput(0)
    myUtil("      The maximum value is " + str(MaxVal) + ".")

    # Gett he NoData value for the Raster.
    NoData = arcpy.Describe(ThisRas).noDataValue
    myUtil("\n     **** NoData value is: " + str(NoData))


    rasterObj = arcpy.Raster(ThisRas)
    NoData2 = rasterObj.noDataValue
    NoDataVal = str(NoData2)
    myUtil("\n     **** NoData2 value is: " + NoDataVal)

    # Check to see if the raster values within the AOI indicate data presence.
    # A value of "0" indicates no data; anything else -> data present.
    if MaxVal == "0":
        myUtil("        Raster %s DOES NOT CONTAIN data within the AOI.\n" % (ThisRas))
    else:
        myUtil("        Raster %s CONTAINS data within the AOI." % (ThisRas))
        outFile = os.path.join(NewFldr,ThisRas)
        arcpy.Copy_management(ThisRas, outFile)
        myUtil("          Raster %s has been copied to saved folder.\n" % (ThisRas))    

myUtil("\nEnd of Script.")
  • Welcome to GIS SE! We're a little different from other sites. We're a Q&A site, not a discussion forum. For questions that involve code we ask that you show us where you are stuck with your own code by including a code snippet in your question. There is an edit button beneath your question which will enable you to do that and a {} button that enables you to format any highlighted code nicely. Please check out our short tour for more about how the site works. Thanks. – PolyGeo Sep 12 at 23:38
  • Do you have a spatial analyst license? Zonal statistics (as table) would fit the bill overlaying your polygons over the raster and returning the maximum. If you do not have spatial analyst then you will need to convert your raster into polygons and perform a union/intersect/identity operation depending on whether you want all your output polygons to be retained or just the ones that intersect your value 1 polygons - but first set the NoData value of your raster to your background value to avoid creating irrelevant polygons. – Michael Stimson Sep 12 at 23:54
  • Thank you PolyGeo - I had forgotten to insert my code and obviously, I'm not very good at inserting it into my message either. Michael - Yes, I do access to SA and was attempting to do this without calling a metered extension if I can help it. – Layne Seely Sep 13 at 16:18
  • Your edit removed the graphics. You probably want to add these back in. – Vince Sep 13 at 16:52
  • Thank you Vince – Layne Seely Sep 16 at 18:25

Your Answer

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

Browse other questions tagged or ask your own question.