1

I have 8,760 (24*365) GRID rasters in a directory. I'd like to use use FocalStatistics analysis to reassign cells to the maximum of a 5x5 window/kernel and extract the resultant values at points specified in a points feature class.

from arcpy import CheckOutExtension
    from arcpy.sa import FocalStatistics, NbrRectangle

    if CheckOutExtension("spatial") != 'CheckedOut':
        raise Exception('Failed to check out spatial analyst extension.)

    # Testing subset of raster data. For simplicity.
    HRRR_rasters = [r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010000',
                    r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010100',
                    r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010200',
                    r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010300',
                    r'E:\HRRR_Data\2014\GRID\1 - JAN_2014\g01010400']

    neighborhood = NbrRectangle(width=5, height=5, units="CELL")
    extract_points = r'C:\Users\jking\Documents\ArcGIS\Default.gdb\sample_points'


    for raster in HRRR_rasters:
        raster_obj = FocalStatistics(raster, neighborhood, 'MAXIMUM', "DATA")

        # ------------------------Psudocode ------------------------------------
        # ExtractValuesToPoints(in_point_features, in_raster, out_point_features) 

        # arcpy.da.InsertCursor to append to master table.
        # Allow raster object to be garbage collected or deleted.

The main problem that I'm having is that ExtractValuesToPoints takes a raster_layer as input rather than a raster object. I have looked into saving the raster to the 'in_memory' workspace without any luck. I've also looked into the properties of the raster object finding that raster_obj.catalogPath & raster_obj.path are None. MakeRasterLayer takes the path and name of the raster dataset which also doesn't work. I can, albeit reluctantly, save the raster to disc, use it for processing and delete it manually. I would like to avoid this write heavy approach though.

  • Try arcpy.Raster("pathToRasterDataset") – FelixIP Apr 16 '15 at 21:49
  • 1
    if you use arcpy.env.overwriteoutput you can overwrite a temporary file and/or arcpy.Delete_management to clean it up when you're finished with it. Another handy approach is to test for existence using if arcpy.Exists(TempRaster): arcpy.Delete_management(TempRaster) on the way in to ensure the raster doesn't exist before creating it. Temp raster objects are GRID data in your %Temp% (or was that %tmp%) directory - they are file system objects nothing 'magic' about them, only that arcpy knows how to handle them. – Michael Stimson Apr 16 '15 at 23:16
2

This usually works:

from arcpy import env
env.overwriteOutput = True
env.workspace = outFolder

raster_obj=FocalStatistics(raster, neighborhood, 'MAXIMUM', "DATA")
raster_obj.save("outZ")
del raster_obj
ExtractValuesToPoints(in_point_features, outFolder+os.sep+"outZ", out_point_features)
  • Add arcpy.Delete_management(outFolder+os.sep+"outZ") to clean up and it's a good answer Felix. – Michael Stimson Apr 16 '15 at 23:18
  • 1
    Yes, outside loop at the very end. Overwrite output option handles it Ok otherwise – FelixIP Apr 17 '15 at 1:00
  • There is only one raster but you're right.. if you're doing many then just before exiting is the best place to put that. – Michael Stimson Apr 17 '15 at 1:03

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.