3

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(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)?

Improve this question
2
  • 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
    Commented Apr 3, 2019 at 7:12
  • Could be an issue with your geometry. You could try the buffer with a zero width trick. Sometimes that cleans up topology. I would also read the data using sf::st_read as it is notably faster and then will let you try @Spacedman recommendation of using an sf class with st_intersection. I would also check if it is multipart geometry. In the sf class, if you print the object the geometry field will tell you if it is POLYGON or MULTIPOLYGON. If multipart I would explode the features before running an intersection.
    – Jeffrey Evans
    Commented Jun 29, 2021 at 18:20

1 Answer 1

Reset to default
1

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

Could you try:

library(cleangeo)
huc = clgeo_Clean(huc)
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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