Clipping/ Zooming in on shapefile using R

Question

I'm working on a shapefile of all the villages of Portugal from this link.
I've uploaded it to R:
#read in shapefile of portugal
portugal = readOGR(dsn = "./communities", layer = "Cont_AAD_CAOP2017", verbose=FALSE)

And I've got the coordinate reference system information
# get coordinate reference system information
st_crs(portugal)

Coordinate Reference System:
  User input: ETRS89 / Portugal TM06 
  wkt:
PROJCRS["ETRS89 / Portugal TM06",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["Portugual TM06",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",39.6682583333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-8.13310833333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Topographic mapping (medium scale)."],
        AREA["Portugal - mainland - onshore."],
        BBOX[36.95,-9.56,42.16,-6.19]],
    ID["EPSG",3763]]

I now want to create a second shapefile using only a region of Portugal, as delimited by some coordinates, namely
W 8°46'00"-- W 7°40'00" / N 41°08'00" -- N 40°39'00"

How can I do this?
Once again, I'm using R

Elio Diaz · Answer

In the following code, first the level 1 administrative boundaries of PRT are downloaded (the provided link did not work) and converted to a sf object; then we convert the character coordinates to dms (degree-minute-second) and to decimal degrees with as.numeric. A polygon is built using those coordinates and afterwards the boundaries object is clipped (intersected) with the bbox polygon:
sf is the modern geospatial library for vector data in R, compatible with the pipe %>% operator; furthermore, it reads vector data way faster than readOGR with read_sf
library(sf)
library(sp)
library(raster)
library(dplyr)

prt = raster::getData(name = "GADM", country = "PRT", level = 1) %>% st_as_sfc()
coords = c("-8°46'00"",
"-7°40'00"",
"41°08'00"",
"40°39'00"")
coords = char2dms(coords, chd = "°") %>%
 as.numeric()

bb_poly = st_bbox(c(xmin = coords[1], xmax = coords[2], ymin = coords[3], ymax = coords[4])) %>% 
  st_as_sfc() %>% st_set_crs(4326)

st_intersection(prt, bb_poly) %>% plot()

Clipping/ Zooming in on shapefile using R

One Answer

Add your own answers!

Ask a Question