How to create a subset of a shp file, with all its properties

I am new to programming in R and with .shp files.

I am trying to take a subsample / subset of a .shp file that is so big, you can download this file from here: https://www.ine.es/ss/Satellite?L=es_ES&c=Page&cid=1259952026632&p=1259952026632&pagename=ProductosYServicios%2FPYSLayout (select the year 2021 and then go ahead).

I have tried several things but none of them work, neither is it worth passing it to sf because it would simply add one more column called geometry with the coordinates listed and that is not enough for me to put it later in the leaflet package.

I have tried this here but it doesn't work for me:

myspdf = readOGR(getwd(), layer = "SECC_CE_20210101") #It works

PV2 = myspdf[myspdf@data$NCA == 'País Vasco', ] #Dont work
PV2 = myspdf[,myspdf@data$NCA == 'País Vasco'] #Dont work

What I intend is to create a sample of myspdf (with data, polygons, plotorder, bbox and proj4string) but I don't want it from all the NCA values (myspdf@data$NCA), I only want those in which data$NCA are 'País Vasco'

In short, I would like to have a sample for each value of the different NCA column.

Is that possible? someone can help me on this? thank you very much.

I have tried this too but the same thing as before appears to me, all 18 variables appear and all are empty:

Pais_V = subset(myspdf, NCA == 'País Vasco')
dim(Pais_V)

Solution

Here's one approach:

library(rgdal)

dlshape=function(shploc, shpfile) {
  temp=tempfile()
  download.file(shploc, temp)
  unzip(temp)
  shp.data <- sapply(".", function(f) {
    fp <- file.path(temp, f)
    return(readOGR(dsn=".",shpfile))
  })
}

setwd("C:/temp")
x = dlshape(shploc="https://www2.census.gov/geo/tiger/GENZ2020/shp/cb_2020_us_aitsn_500k.zip", "cb_2020_us_aitsn_500k")

x<-x$.  # extract the shapefile

mycats<-c("00","T2","T3","28")
x2<-subset(x,  x$LSAD %in% mycats)  # subset using the list `mycats`
mypal=colorFactor("Dark2",domain=x2$LSAD)

library(leaflet)
leaflet(x2) %>% addPolygons(weight=.2, color=mypal(x2$LSAD))

dlshape function courtesy of @yokota