I am trying to do something relatively simple but with large shapefiles. I have a shapefile of protected areas in the USA:
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(pas)
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
Intersection:
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)?
sf
and see what happens. – Spacedman Apr 3 '19 at 7:12