I am trying to do something relatively simple but with large shapefiles. I have a shapefile of protected areas in the USA: protected areas

and watershed boundaries: Watershed boundaries

I want to calculate the percentage and proportion of protected areas in each watershed. I want to retain the designation of each protected area as well (i.e., federal, state, etc.)

To do this I started to try to intersect the two shapefiles:

pas <- readOGR(dsn = ".", layer = "PADUS1_4Combined") # protected areas

huc <- readOGR(dsn = ".", layer = "HUC2_USA_PRJ") # hydrologic units (2 digit, USGS)

Checked the CRS



crs(huc) CRS arguments: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

crs(pas) CRS arguments: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0


test <- gIntersection(pas, huc)

I get this error:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false, : rgeos_crdMat2LinearRing: linearRing not created

I searched for this error and cannot find anything. I am dealing with a very large dataset here so I am wondering if this is the issue.

Is there a workaround using simple features (sf)?

  • It could be the size of the dataset (you've not said how big it is) or it could be a problem with the geometry of a feature - perhaps a polygon not valid or something. Try it with subsets. Also try sf and see what happens. – Spacedman Apr 3 at 7:12

It may be there is a problem with your shapefiles which could be cleaned before you attempt the intersection?

Could you try:

huc = clgeo_Clean(huc)

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.