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.

Improve this question
2
  • Try arcpy.Raster("pathToRasterDataset")
    – FelixIP
    Apr 16, 2015 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, 2015 at 23:16

1 Answer 1

Reset to default
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)
Improve this answer
3
  • Add arcpy.Delete_management(outFolder+os.sep+"outZ") to clean up and it's a good answer Felix.
    – Michael Stimson
    Apr 16, 2015 at 23:18
  • 1
    Yes, outside loop at the very end. Overwrite output option handles it Ok otherwise
    – FelixIP
    Apr 17, 2015 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, 2015 at 1:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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