1

I am having trouble with R's raster::intersect method. I give intersect a SpatialPolygonsDataFrame and SpatialPoints as documented in the raster reference manual (p. 111). The result acts like it ignores second parameter. My code follows.

Make spatial points, note pcts_10_sp is a SpatialPolygonsDataFrame

temp_coords <- coordinates(pcts_10_sp) 
temp_crs <- crs(pcts_10_sp)
sp2010 <-SpatialPoints(temp_coords , proj4string=temp_crs )
summary(sp2010)

Results of SpatialPointsDataFrame summary

Object of class SpatialPoints
Coordinates:
                min       max
coords.x1 -97.03242 -96.52624
coords.x2  32.55765  32.98299
Is projected: FALSE 
proj4string : [+proj=longlat +datum=NAD27+ellps=clrk66+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat]
Number of points: 694

Read polygon shapefile and project

pcts_12_sp <- rgdal::readOGR(dsn="Precinct12.shp", layer="Precinct12", verbose = TRUE) 
pcts_12_sp <- sp::spTransform(pcts_12_sp, CRS("+proj=longlat +datum=NAD27"))
summary(pcts_12_sp)

Results of SpatialPolygonsDataFrame summary

Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x -97.03810 -96.5166
y  32.54508  32.9892
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=NAD27 +ellps=clrk66
+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat]
Data attributes:
    2012_PCT         LD     
 1000   :  1   100    : 72  
 1001   :  1   114    : 68  
 1002   :  1   103    : 66  
 1003   :  1   108    : 66  
 1004   :  1   109    : 61  
 1005   :  1   105    : 56  
 (Other):768   (Other):385 

Run raster:intersect function

RasterShp <- raster::intersect(pcts_12_sp, sp2010)
summary(RasterShp)

Results of intersect summary

Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x -97.03810 -96.5166
y  32.54508  32.9892
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=NAD27 +ellps=clrk66
+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat]
Data attributes:
    2012_PCT         LD     
 1000   :  1   114    : 57  
 1001   :  1   100    : 55  
 1003   :  1   109    : 52  
 1004   :  1   108    : 49  
 1005   :  1   103    : 46  
 1006   :  1   111    : 45  
 (Other):610   (Other):312  

The barebones plot of polygons without any points.

plot(RasterShp) 

Plot result

  • What are you expecting and where do you think the function is going wrong. Since we do not have access to your data (polygon) it is difficult to troubleshoot. In this case of intersect your points would receive the polygon attributes but I do not think that the polygons would be effected. For this you would have do use sp::over. – Jeffrey Evans May 27 '16 at 21:06
  • If you raster::intersect polygons with points it will actually intersect the polygon with the extent of the points, i.e. intersect(poly, pts) is the same as intersect(poly, extent(pts)). Looking at the extent of your objects, it appears the points have a slightly smaller extent than the polygons, so I would expect the polygons to be trimmed at the edges. I suspect you're trying to get the subset of polygons that contain a point, in which case you want overlays from the sp package. There's a great vignette on this. – Matt SM May 28 '16 at 20:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.