Updating raster value to reflect total number of points within cell using Python?

Question

I have to do a task in Python/ArcPy/GDAL/whatever because I'm trying to automate this process for a hundred or so different raster/shapefile combos.

This is what I need to do:

I have an empty raster and a shapefile with thousands of points. I would like to update the value of each cell/pixel in the raster to reflect the number of times a unique point lands within that cell.

At first I was thinking something like this:

for cell in raster:
    for point in shapefile:
         if point intersects with (/is contained in) cell: cell_value + = 1

However, I don't know how to implement this because I don't see how to intersect a raster and a polygon. Converting the raster into polygons is out of the question because the rasters are very large with a lot of cells.

I'm not married to the idea of looping through both the individual cells and points because I feel like it's unnecessary, so if there is an ArcGIS tool(s) that I could manipulate, and/or another method in python I could use to reach my goal, I'm more than happy to try that out.

Answer 1

Using any of the tools you mentioned you can preprocess the points so that one feature (with a count) occurs in each raster cell. Converting that to a raster using the standard built-in method completes the task.

The advantage of this method is that it provides complete control over the conversion process, so you know exactly what is happening. Moreover, it easily generalizes: with only trivial modifications (at the summarize step) you can produce rasters whose values are any statistical summary of the points falling within each cell (of which the count is a simple example).

The preprocessing is relatively simple: you need to identify the cell in which any point (x,y) falls. The raster's extent is given by an origin (Ox, Oy) and corresponding cell sizes (h, k) (usually h=k). The (column, row) coordinates of (x,y)'s cell are

(i,j) = floor( (x-Ox)/h, (y-Oy)/k )

That requires two quick calculations to create new shapefile attributes i and j. Summarize the attribute table on the concatenation of i and j (which will identify the cells), retaining the total count (or any other summary) for each unique (i,j) pair, along with the values of i and j themselves. Convert the (i,j) attributes in the summary table back to coordinates of cell centers as

(x,y) = ( Ox + h*(i+1/2), Oy + k*(j+1/2) )

This produces a table with three attributes: x, y, and the summary value for the points corresponding to (x,y)'s cell. Proceed in any way you prefer to convert this into a raster.

Answer 2

I suggest using the Point to Raster tool. With default settings, it chooses one point per cell to generate raster values. Using the cell_assignment option of COUNT, however, the raster stores how many points occurred in each cell.

This produces new rasters rather than working with existing ones, but the process of combining rasters for analysis is simpler than combining raster/vector data.

Example arcpy code (see help link above for full context) would be:

# Set local variables
inFeatures = "ca_ozone_pts.shp"
valField = "ELEVATION"
outRaster = "c:/output/ca_elev02"
assignmentType = "COUNT"
priorityField = ""
cellSize = 2000

# Execute PointToRaster
arcpy.PointToRaster_conversion(inFeatures, valField, outRaster, 
                               assignmentType, priorityField, cellSize)