Fastest way to test if raster contains data within a polygonal area

Question

I am writing a python script for a geoprocessing tool in arcmap. I have a polygon in WGS 1984 Web Mercator (Auxiliary Sphere) and a raster grid, raster1, in WGS 1984. I reproject the polygon to create a mesh of points that are contained within the boundary of the polygon and laid out in a uniform grid. Now, the goal is to check raster1 at each point in the mesh and detrimine if the corrosponding cell has a value at that xy. If it does remove the point from the mesh. If it doesn't, keep the point. In the end I want to have a fragment of the mesh of points which contains all points for which the raster had no value.

My first attempt was to loop through the points and use arcpy.GetCellValue_management() to see what the value was at each point. But the GetCellValue_management fuction seems to take a painfully slow 2-3 seconds (I have a lot of points). Anyone know why this takes so long? An identify on the same raster in arcmap is almost instant!

So now I'm thinking about using Clip_management to clip out the intersection of the polygon and the raster and then RasterToPolygon to get a polygon representing that intersection. Then finally using Erase_analysis to erase from the point mesh those points in the intersection.

What do you guys think? Is plan two good? Whats the fastest way of doing this?

Answer 1

The fastest and easiest solutions usually exploit simple, fast, built-in procedures in the GIS.

If I have understood correctly, you begin with a polygon and a source raster of data (both shown projected):

Values of the source raster are shown in shades of red. Many of the source cells are missing data (the yellow background shows through).

Exploit the fast rasterization capabilities to create a raster "fishnet" of the polygon. Here, I have specified a much coarser cellsize than used by the source grid (and I have offset the origin, too, so that there is no nice relationship between this new grid and the source grid):

(The shading doesn't mean anything: it is varied randomly as a quick way to show the cells.)

Now you merely have to convert the source raster to another raster using the origin and cellsize of the rasterized polygon. This requires resampling. Specify nearest neighbor resampling. (This is the default in ArcGIS I believe.) To the center of each coarse polygonal cell, this assigns the value of the nearest (i.e., overlaying) source grid. In so doing, mask the resampling to the polygon. This rapid built-in operation does all the work:

The green cells are those for which source data are available.

Finally, a little bit of map algebra polishes the result. In pseudocode (which I'll use because syntax varies among implementations), the operation is

SetNull(Not(IsNull([Source])), 1).

This places indicator values (1's, or anything you like) at all coarse cells receiving NoData values from the source grid. (You can probably skip the previous step and do this directly, because the resampling will happen automatically.)

Here you can see the source grid overlaid by the result, again shown in grayscale to distinguish individual cells. If needed, convert this result to a layer of vector points (one point per cell center). Of all the possible raster-vector conversions you could do in this analysis, this will be one of the fastest.

Let's recapitulate: the workflow requires

1. Stipulating the origin and cellsize of the desired output grid.
2. Stipulating that the polygon will be used to mask raster operations.
3. Resampling the source grid (which may be optional, depending on the GIS).
4. One fast map algebra command.
5. Optionally, conversion of the resulting raster to vector format.

You will likely find the overall computation time to be faster than creating a vector "fishnet" layer in the first place, which would only get you through step 2 anyway.

Answer 2

My method is a bit roundabout, and probably wrong in at least one sense, but if it works for you, go for it.

Basically, I'd do the following:

1. Convert the raster to a polygon with the value in a table column - seems like both our methods use this step.
2. Spatial join the newly created polygon to the points so that one layer has all of the information.
3. Simply use the cursor over the points layer and the data is already transferred. In my experience, a spatial join for a large number of features to a large number of features and then querying values and utilizing them is pretty quick. You should be able to use an UpdateCursor to remove the points quickly, then just save it at the end.

I hope that helps you - sounds like your method will work too - not sure what will be faster.

Answer 3

I think it might be quicker to use rasters for all of the calculations. For example, you could create a raster from your polygon. Then, identify all of the areas that touch the raster. You can do this through reclassification, and addition of values. You could set the polygon raster to be 0, the orginial raster to 1 and add them together. Then any cell that is 0 (after this process) is an area that does not touch the original raster.