7

I have two datasets: point feature class and raster (stored in file geodatabase).

The spatial extent of the raster is smaller than the point dataset.

Is there a way I could select (or create new, temporary layer/feature class) only the points that fall into the extent of the raster?

Basically I'd like to have an output similar to Clip or Intersect tools, but they don't seem to accept rasters as inputs.

Is there readily available arcpy tool(s) that could be used to achieve that? Or would it require some custom solution?

2

4 Answers 4

6

The attached script which utilizes Raster Domain will work for you.

# Import arcpy module
import arcpy
from arcpy import env

# Check out any necessary licenses
arcpy.CheckOutExtension("3D")

# Set Workspace
env.workspace = r"C:\data"

# Create a polygon around the raster boundary
arcpy.RasterDomain_3d("inRaster.tif", "outShp", "POLYGON")

# Intersect points with boundary polygon
arcpy.Intersect_analysis(["outShp.shp", "points.shp"], "intersect_output", "ALL", "", "INPUT")

For a direct, non ESRI approach, check out Geospatial Modelling Environment's (GME) isectpntrst (Intersect Points With Raster).

0
6

You could use the tool "Raster to Polygon" from the "Conversion Tools". As output you get a shapefile with the extent of your raster image. With this shapefile and the other one you could run the intersect command.

import arcpy
from arcpy import env
env.workspace = "C:/data"
arcpy.RasterToPolygon_conversion("zone", "c:/output/zones.shp", "NO_SIMPLIFY",
                                  "VALUE")
1
  • Thanks. Was thinking about this option as well and it indeed might be a viable solution.
    – radek
    Aug 28, 2012 at 13:29
5

Use these tools from ArcToolbox:

1-Spatial Analyst Tools> Raster Creation> Create Constant Raster ("Constant value" can
  be 0.5, 1, ... and "Output extent" must be your overlay raster )

2-Conversion Tools> From Raster> Raster to Polygon (Use output raster from step1)

3-Analysis Tools> Overlay> Intersect (the polygon with the points)

4-Done!

and here is the model:enter image description here

1

Raster to Polygon > Intersect will certainly work.

Depending on the size of your raster (in terms of rows/columns, not geographic extent) the Raster to Polygon step can be computationally expensive.

If your raster doesn't contain any areas of nodata, or you don't care if you include points falling on nodata, a very fast solution would be to set the extent environment to the extent of the raster, and then create a copy of your file which will only contain points within the extent of your raster's bounding box. The 'Copy' tool doesn't use the extent environment, so I'd use the 'merge' tool with only the point dataset as an input.

import arcpy
from arcpy import env

#set environment to raster extent
env.extent = yourraster

arcpy.Merge_management(yourpoints, output_points)

The big value to this method is that it is very fast.

If your raster contains areas of nodata that you want to exclude from your final point dataset, you can use the Spatial Analyst tool "Extract Values to Points". Again, set your extent environment variable to the extent of the raster, and then use the ExtractValuesToPoints tool to create a point file of those points that fall within the extent of your raster with a new field called RASTERVALU. Points falling on nodata bins will have a value of '-9999' (warning: if your raster contains real values of -9999, that could be an issue - not sure if there's a workaround).

import arcpy
from arcpy import env
from arcpy import sa

#Check out Spatial Analyst license or raise exception
if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.CheckOutExtension("Spatial")
else:
    sys.exit('Spatial Analyst license not available.')

env.extent = yourraster

sa.ExtractValuesToPoints(yourpoints, yourraster, outputpoints)

Then you can select out those points with nodata values to remove them from your output point file, or extract the relevant points via whatever method works for you.

Bonus to this is that it attaches the values of your raster to each of your points, so you can not only select the points that don't fall on nodata, but you can display and manipulate your other points relative to the associated raster values. NOTE: ExtractValuesToPoints does not work with multipart point features.

ExtractMultiValuesToPoints would also work.

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.