I have an issue. I have a set of points (n=~2000) and I would like to extract the pixel values from multiple rasters (a single-raster solution is also ok) for location of that point as well as that of all eight neighbours.
Here is how I would consider doing this if I absolutely had to have every neighborhood value (which is what it sounds like). First create 2 fields for the point geometry (one for x and one for y) and calculate the geometry for the x field and y field. I would be sure my point file was in a projected coordinate system like UTM Meters or State Plane Feet. Better if the units match units represented in the cell size of the raster. Then I'd create 16 more fields ULX, ULY, UCX, UCY, URX, URY, CLX, CLY, CRX, CRY, LLX, LLY, LCX, LCY, LRX, LRY.
I'm going to describe my example using 30 meter cells. So next I'd use the x and y fields for the point to create new x and y values for ULX, LCX, and LLX calculate each of those fields with the field calculator X - 30. for UC and LC it is just the same value as X, and for UR, CR, and LR it is x + 30. Then to calculate the Y values CL and CR are equal to Y. The top row UL, UC, and UR are Y + 30 and LL, LC, LR are Y - 30. So now all the fields are populated. I'd create 8 new point feature classes or shapefiles using the 16 x and y fields using "Add XY Data". Keeping an ID value the same in all of them (FID), you can then run "Extract Values to Points" or "Sample" from the Spatial Analyst Toolbox on the original file (for the central cell) and all the other point feature classes. These results could then be brought individually into a spreadsheet and put back together so that each original ID has a value for UL, UC, UR, CL, Original Point, CR, LL, LC, and LR.
You could script the adding of the fields and calculating of the values if you needed to do this again for different sets of points.
For a 5 x 5 neighborhood, you would need to be more creative in naming your fields and there would be 25 total feature classes / shapefiles, but it could be done this way too.
You appear to be looking at a suite of software options. A way of doing it in ArcMap model builder using off the shelf tools could be:
- point to raster (ensure snap to raster environment is set.)
- Expand (by 1 pixel to create your block of nine)
- Extract by mask.
This method assumes that your points are not so close that their masks overlap. If overlap is going to be an issue then you could use a feature iterator to feed out single points and process them one at a time, this would generate a raster per point.