4

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.

  • If you convert points to raster, you can do it in a way that counts the points in that cell (by using COUNT for the cell_assignment option). Would it be feasible for you to convert all your point shapefiles into rasters, then analyze the raster values (which will be point counts)? – Erica Apr 4 '14 at 15:19
  • @Erica thanks for responding! I'm not sure if that would be feasible because I've never used the tool. I'm in the process of trying to produce my shapefiles, but I'm expecting the number of points to be as high as 70 million. What do you think? – user20408 Apr 4 '14 at 15:24
  • I would give it a shot; it can be used in arcpy (ref. the link in my prev. comment for documentation and an example, just use assignmentType = "COUNT" instead of "MAXIMUM"). I honestly don't know how its speed would compare to counting. It will be making a raster for each of your point shapefiles, so storage may become a concern (although you can delete the raster after you're done with it). – Erica Apr 4 '14 at 15:37
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)
  • I used this tool in the way you said and it seemed to have worked. Except for the fact that there are NoData values present in large 'chunks' around the points, in addition to the correct values of 1's and 2's where the points were. I suppose I can add another step to fill in the no data values with 0's, but do you know why this is happening, or what I could do to fix it? Thanks! – user20408 Apr 4 '14 at 16:57
  • 1
    You can use the CON tool to set the NOATA values to a zero whilst keeping all other cells the same value. – Hornbydd Apr 4 '14 at 23:13
  • 2
    It's presumably doing that because the "default" point-to-raster process assigns values, not counts, and so NoData would be completely appropriate in that situation. The fix is quick with Con, like @Hornbydd said -- there is an example to convert NoData to 0 in the Con help page (example 2). You might not need to do this, it depends on what subsequent analysis you want to do with the rasters; however, it's straightforward to script and automate if you need the 0's. – Erica Apr 5 '14 at 0:31
3

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.

  • Thanks whuber. This looks to be the route I'm going to have to take, as I can no longer use arcpy or arcgis, which means I can't use the point to raster tool. So now, I have a csv file of lat, long coordinates, and I need to produce a raster at 3 varying resolutions (25km, 5km, and 1km) that depicts how many points fall into each cell. I have to do this in gdal now. Do you have any suggestions for implementing this algorithm using gdal? Thanks! – user20408 Dec 2 '14 at 16:39
  • @user20408 I would do this in R and take advantage of packages that provide GDAL access, such as rgdal. It would take about a dozen lines of code--one to read the file, a few to perform the calculations, and a few more to create a raster and save it as a dataset. Alternatively, if you have Python-GDAL bindings (I'm not familiar with this so I can't give details), you could do all the work directly in Python. – whuber Dec 2 '14 at 17:12

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.