Exercises for the meeting of May 10, 2019
Violeta Sosa
May 09, 2019
Exercises for the meeting of May 10, 2019
Exercises related to simple features in R (1). Try to upload your answer to the LearnWeb assigned as a .zip file (not a .rar file) containing two files: (i) an R markdown file and (ii) an html file that resulted from knitting the R markdown file. (If this is problematic, use pdf or something else.)
##Exercise 1 Look for a shapefile; download it; import it with
sf::st_read, and plot it. What does it show?
shapeName <- st_read('depto/depto.shp')## Reading layer `depto' from data source
## `C:\Users\viole\Documents\SDR_geotech\Assigment 4\depto\depto.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 165298.4 ymin: 20565.69 xmax: 1804280 ymax: 1984871
## Projected CRS: Bogota 1975 / Colombia Bogota zoneshapeName## Simple feature collection with 33 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 165298.4 ymin: 20565.69 xmax: 1804280 ymax: 1984871
## Projected CRS: Bogota 1975 / Colombia Bogota zone
## First 10 features:
## DPTO NOMBRE_DPT AREA PERIMETER HECTARES EFPROM
## 1 5 ANTIOQUIA 63351855547 1963728.8 6335185.6 0.9668
## 2 8 ATLANTICO 3360765350 240936.2 336076.5 0.8347
## 3 11 SANTAFE DE BOGOTA D.C 1650947779 323322.5 165094.8 1.0000
## 4 13 BOLIVAR 26141894528 1309428.0 2614189.5 0.6016
## 5 15 BOYACA 23352582464 1364539.9 2335258.2 0.5466
## 6 17 CALDAS 7558199876 603282.5 755820.0 0.7112
## 7 18 CAQUETA 90180868829 1888506.9 9018086.9 0.6604
## 8 19 CAUCA 29742787301 1243389.0 2974278.7 0.5588
## 9 20 CESAR 22973095680 1080343.7 2297309.6 0.5219
## 10 23 CORDOBA 25059485606 814093.4 2505948.6 0.9833
## geometry
## 1 MULTIPOLYGON (((754872.1 14...
## 2 MULTIPOLYGON (((913511.2 16...
## 3 MULTIPOLYGON (((1006436 102...
## 4 MULTIPOLYGON (((881889.9 16...
## 5 MULTIPOLYGON (((1206428 126...
## 6 MULTIPOLYGON (((931978.6 11...
## 7 MULTIPOLYGON (((931968.4 76...
## 8 MULTIPOLYGON (((540370.4 82...
## 9 MULTIPOLYGON (((1087325 169...
## 10 MULTIPOLYGON (((809020.2 15...plot(shapeName, main="Colombia shapefile")
The plot shows a map per each one of the attributes defined in the shapefile’s attribute table.
##Exercise 2 Try plotting one of the attribute variables using a different color scheme
shapeName%>%select(AREA) %>%plot(graticule = TRUE, axes = TRUE)
shapeName%>%select(AREA) %>%plot(graticule = TRUE, axes = TRUE, col = heat.colors(10))
shapeName%>%select(DPTO) %>%plot(graticule = TRUE, axes = TRUE, col = topo.colors(50))
##Exercise 3 What is the class of the object plotted? Plot only the
geometry using st geometry
class(shapeName)## [1] "sf" "data.frame"plot(st_geometry(shapeName), main="Colombia Geometry")
##Exercise 4 Which layers does the GeoPackage http://www.geopackage.org/data/sample1_2.gpkg have?
Import them both in R, using sf. Call the object with
feature geometries x1.
st_layers("sample1_2.gpkg")## Driver: GPKG
## Available layers:
## layer_name geometry_type features fields
## 1 counties Multi Polygon 3141 9
## 2 countiestbl NA 3141 11dataset <- ("sample1_2.gpkg")
counties<-read_sf(dataset, "counties")
countiestbl<-read_sf(dataset, "countiestbl")
x1<-counties
x1## Simple feature collection with 3141 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -178.2176 ymin: 18.92179 xmax: -66.96927 ymax: 71.40624
## Geodetic CRS: NAD83
## # A tibble: 3,141 x 10
## NAME STATE_NAME STATE_FIPS CNTY_FIPS FIPS AREA POP1990 POP2000 POP90_SQMI
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <int>
## 1 Lake ~ Minnesota 27 077 27077 1784. 4076 4651 2
## 2 Ferry Washington 53 019 53019 2280. 6295 7199 3
## 3 Steve~ Washington 53 065 53065 2530. 30948 40652 12
## 4 Okano~ Washington 53 047 53047 5306. 33350 38640 6
## 5 Pend ~ Washington 53 051 53051 1445. 8915 11752 6
## 6 Bound~ Idaho 16 021 16021 1279. 8332 10068 7
## 7 Linco~ Montana 30 053 30053 3746. 17481 18859 5
## 8 Flath~ Montana 30 029 30029 5232. 59218 73438 11
## 9 Glaci~ Montana 30 035 30035 3124. 12121 12626 4
## 10 Toole Montana 30 101 30101 1943. 5046 4556 3
## # ... with 3,131 more rows, and 1 more variable: Shape <MULTIPOLYGON [°]>x1%>%select(AREA) %>%plot(graticule = TRUE, axes = TRUE, main="Colombia (Only) Geometry Plot")
##Exercise 5 Create an sf object with a point geometry
at latitude 36.21 and longitude -81.19
(sfc = st_sfc(st_point(c(-81.19,36.21)), crs = 4269))## Geometry set for 1 feature
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -81.19 ymin: 36.21 xmax: -81.19 ymax: 36.21
## Geodetic CRS: NAD83## POINT (-81.19 36.21)##Exercise 6 Query the features of x1 at this point: plot the resulting geometry, and show the attributes associated with this geometry.
geometry=st_intersects(sfc,x1)
geometry## Sparse geometry binary predicate list of length 1, where the predicate
## was `intersects'
## 1: 1874plot(st_geometry(x1[unlist(geometry),]),expandBB = c(2,2,2,2), col=rgb(red = 1, green = 0, blue = 0, alpha = 0.5), main="Resulting Geometries")
plot(st_geometry(x1), add=TRUE)
plot(sfc, expandBB = c(10, 10, 10, 10), col = "red", lwd = 5, add=TRUE)
x1[unlist(geometry),]## Simple feature collection with 1 feature and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.54593 ymin: 35.98946 xmax: -80.87066 ymax: 36.43128
## Geodetic CRS: NAD83
## # A tibble: 1 x 10
## NAME STATE_NAME STATE_FIPS CNTY_FIPS FIPS AREA POP1990 POP2000 POP90_SQMI
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <int>
## 1 Wilkes North Carolina 37 193 37193 767. 59393 64307 77
## # ... with 1 more variable: Shape <MULTIPOLYGON [°]>##Exercise 7 For the nc dataset that ships with package
sf, compute the area of county Columbus
(file <- system.file("gpkg/nc.gpkg", package="sf"))## [1] "C:/Users/viole/OneDrive/Documentos/R/win-library/4.1/sf/gpkg/nc.gpkg"(file %>% read_sf() -> nc)## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
## # A tibble: 100 x 15
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
## <dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.114 1.44 1825 1825 Ashe 37009 37009 5 1091 1 10
## 2 0.061 1.23 1827 1827 Alle~ 37005 37005 3 487 0 10
## 3 0.143 1.63 1828 1828 Surry 37171 37171 86 3188 5 208
## 4 0.07 2.97 1831 1831 Curr~ 37053 37053 27 508 1 123
## 5 0.153 2.21 1832 1832 Nort~ 37131 37131 66 1421 9 1066
## 6 0.097 1.67 1833 1833 Hert~ 37091 37091 46 1452 7 954
## 7 0.062 1.55 1834 1834 Camd~ 37029 37029 15 286 0 115
## 8 0.091 1.28 1835 1835 Gates 37073 37073 37 420 0 254
## 9 0.118 1.42 1836 1836 Warr~ 37185 37185 93 968 4 748
## 10 0.124 1.43 1837 1837 Stok~ 37169 37169 85 1612 1 160
## # ... with 90 more rows, and 4 more variables: BIR79 <dbl>, SID79 <dbl>,
## # NWBIR79 <dbl>, geom <MULTIPOLYGON [°]>nc%>%select(NAME)%>%filter(NAME=="Columbus") %>%plot(graticule = TRUE, axes = TRUE, col=rgb(red = 1, green = 0, blue = 0, alpha = 0.5), ADD=TRUE, main="Columbus County")
columbus = nc %>% filter(NAME == "Columbus")
st_area(columbus)## 2450830549 [m^2]##Exercise 8 Compute also the same area after transforming the
geometry to EPSG 2264, and compare the value with that
obtained from the unprojected data; express the difference in a
percentage.
nc_newgeo<-nc
(nc_newgeo %>% st_transform(2264) -> nc_newgeo.2264)## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 406265 ymin: 48359.7 xmax: 3052877 ymax: 1044143
## Projected CRS: NAD83 / North Carolina (ftUS)
## # A tibble: 100 x 15
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
## * <dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.114 1.44 1825 1825 Ashe 37009 37009 5 1091 1 10
## 2 0.061 1.23 1827 1827 Alle~ 37005 37005 3 487 0 10
## 3 0.143 1.63 1828 1828 Surry 37171 37171 86 3188 5 208
## 4 0.07 2.97 1831 1831 Curr~ 37053 37053 27 508 1 123
## 5 0.153 2.21 1832 1832 Nort~ 37131 37131 66 1421 9 1066
## 6 0.097 1.67 1833 1833 Hert~ 37091 37091 46 1452 7 954
## 7 0.062 1.55 1834 1834 Camd~ 37029 37029 15 286 0 115
## 8 0.091 1.28 1835 1835 Gates 37073 37073 37 420 0 254
## 9 0.118 1.42 1836 1836 Warr~ 37185 37185 93 968 4 748
## 10 0.124 1.43 1837 1837 Stok~ 37169 37169 85 1612 1 160
## # ... with 90 more rows, and 4 more variables: BIR79 <dbl>, SID79 <dbl>,
## # NWBIR79 <dbl>, geom <MULTIPOLYGON [US_survey_foot]>columbus_new = nc_newgeo.2264 %>% filter(NAME == "Columbus")
(area4267=st_geod_area(columbus))## 2450241815 [m^2](area2264=(st_area(columbus_new) %>% set_units(m^2)))## 2450348999 [m^2]#New area in percentage
((area2264/area4267)*100)## 100.0044 [1]##Exercise 9 Create a plot with the nc states outlines
that has all states that are partly within 50 km from Columbus filled
with color green.
columbusarea = nc[nc$NAME == "Columbus", ]
(selection= st_is_within_distance(columbusarea,nc, dist = 50))## Sparse geometry binary predicate list of length 1, where the predicate
## was `is_within_distance'
## 1: 94, 96, 97, 98, 100plot(st_geometry(nc), main="States partly within 50km from Columbus County")
plot(st_geometry(nc[unlist(selection),]), add = TRUE, col = 'green')
##Exercise 10 Create a plot with the nc states outlines
that has all states that are entirely within 100 km from Columbus filled
with color red.
(columbusarea = nc[nc$NAME == "Columbus", ])## Simple feature collection with 1 feature and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -79.0745 ymin: 33.94867 xmax: -78.15479 ymax: 34.4772
## Geodetic CRS: NAD27
## # A tibble: 1 x 15
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
## <dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.24 2.37 2232 2232 Columbus 37047 37047 24 3350 15 1431
## # ... with 4 more variables: BIR79 <dbl>, SID79 <dbl>, NWBIR79 <dbl>,
## # geom <MULTIPOLYGON [°]>plot(st_geometry(nc), main="States entirely within 100km from Columbus County")
#Columbus transform for buffer distance calculation
nc_tr = st_transform(columbusarea, 3857)
plot(st_transform(st_buffer(nc_tr, 100000),4267),col=rgb(red = 1, green = 0, blue = 0, alpha = 0.5), border = rgb(red = 1, green = 0, blue = 0, alpha = 0.5), add = TRUE)
otro=st_transform(st_buffer(nc_tr, 100000),4267)
contain_true=st_contains(otro,nc)
plot(st_geometry(nc[unlist(contain_true),]), add = TRUE, col = 'red')