Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this question
Related:… – Chad Cooper Aug 28 '12 at 14:00
@ChadCooper Thanks. Your answer to this question might be indeed good solution. Will give it a go. – radek Aug 28 '12 at 14:19
up vote 6 down vote accepted

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

# 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).

share|improve this answer

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",
share|improve this answer
Thanks. Was thinking about this option as well and it indeed might be a viable solution. – radek Aug 28 '12 at 13:29

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)


and here is the model:enter image description here

share|improve this answer

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":
    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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.