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

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

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)