Warm Up - sf Tutorial

In this document I will show some basic functions from the sf package so you can have the intuition of how it works and what kind of operations can be performed. You will learn how to read geographic data into your R workflow and how to write geographic data in you computer. In addition you will see a workflow analysis and the map making process.
You cand download all files used from my Github Repository. Please notice that all datasets are in the folder /Notebooks/Montevideo_Data.

library(sf)

package ‘sf’ was built under R version 3.5.2Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3

library(tmap)

package ‘tmap’ was built under R version 3.5.2

library(tidyverse)

── Attaching packages ────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.1.0     ✔ purrr   0.2.5
✔ tibble  2.1.1     ✔ dplyr   0.8.1
✔ tidyr   0.8.2     ✔ stringr 1.3.1
✔ readr   1.1.1     ✔ forcats 0.3.0
package ‘tibble’ was built under R version 3.5.2── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

setwd("~/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data")

The working directory was changed to /Users/sabogal/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

#The path to the working directory should change according to where you allocated the repository in your computer. 
#setwd('..')

getwd()

[1] "/Users/sabogal/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks"

---

1 Points (Accidents)

The most common type for spatial datasets are ESRI shape files. You can use the st_read() function and set the dsn (data source name) argument. Conversely, to write an spatial object the st_write() function is defined.

setwd("~/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/Accidentes2006-2010")

The working directory was changed to /Users/sabogal/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/Accidentes2006-2010 inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

Accidents <- st_read("accidentes2006-2010.shp")

Reading layer `accidentes2006-2010' from data source `/Users/sabogal/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/Accidentes2006-2010/accidentes2006-2010.shp' using driver `ESRI Shapefile'
Simple feature collection with 41121 features and 7 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 558003.5 ymin: 6134508 xmax: 588206.1 ymax: 6157857
epsg (SRID):    NA
proj4string:    NA

The above function created the object Accidents. As usual, you can print the object in your console (or rmarkdown document) by typing the name of the object or get a general overview with the str() and summary() functions. An easy way to make a map is by passing the sf object to the plot() function.

Accidents

Simple feature collection with 41121 features and 7 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 558003.5 ymin: 6134508 xmax: 588206.1 ymax: 6157857
epsg (SRID):    NA
proj4string:    NA
First 10 features:
   CALLE CRUCE  TIPO ANIO         NOM_CALLE              NOM_CRUCE
1   4274  3075 FATAL 2010 DR JOSE MARTIRENE CNO GRAL LEANDRO GOMEZ
2   2775  2427 GRAVE 2010 CNO MANUEL FORTET         CNO JOSE DURAN
3   2775  2427 GRAVE 2010 CNO MANUEL FORTET         CNO JOSE DURAN
4   3918  7425  LEVE 2010         AV LEZICA                 YEGROS
5   3918  7425  LEVE 2010         AV LEZICA                 YEGROS
6    126  6213  LEVE 2010      AV AGRACIADA        MARIANO SAGASTA
7    126  6213  LEVE 2010      AV AGRACIADA        MARIANO SAGASTA
8   3162   762  LEVE 2010      MATTO GROSSO   CALDERON DE LA BARCA
9   5115  4794 GRAVE 2010          PARAGUAY             NUEVA YORK
10  3618  4494  LEVE 2010           JOANICO          CIPRIANO MIRO
                                      ESQUINA                 geometry
1  DR JOSE MARTIRENE Y CNO GRAL LEANDRO GOMEZ POINT (575552.3 6145117)
2          CNO MANUEL FORTET Y CNO JOSE DURAN   POINT (572329 6147839)
3          CNO MANUEL FORTET Y CNO JOSE DURAN   POINT (572329 6147839)
4                          AV LEZICA Y YEGROS POINT (570472.3 6148494)
5                          AV LEZICA Y YEGROS POINT (570472.3 6148494)
6              AV AGRACIADA Y MARIANO SAGASTA POINT (571388.7 6142519)
7              AV AGRACIADA Y MARIANO SAGASTA POINT (571388.7 6142519)
8         MATTO GROSSO Y CALDERON DE LA BARCA POINT (570552.1 6147980)
9                       PARAGUAY Y NUEVA YORK POINT (573695.4 6138081)
10                    JOANICO Y CIPRIANO MIRO POINT (578559.9 6140612)

Question: Why are not the 41121?

str(Accidents)

Classes ‘sf’ and 'data.frame':  41121 obs. of  8 variables:
 $ CALLE    : int  4274 2775 2775 3918 3918 126 126 3162 5115 3618 ...
 $ CRUCE    : int  3075 2427 2427 7425 7425 6213 6213 762 4794 4494 ...
 $ TIPO     : Factor w/ 3 levels "FATAL","GRAVE",..: 1 2 2 3 3 3 3 3 2 3 ...
 $ ANIO     : int  2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
 $ NOM_CALLE: Factor w/ 1680 levels "1 DE MARZO","12 DE DICIEMBRE",..: 619 472 472 191 191 134 134 1226 1321 976 ...
 $ NOM_CRUCE: Factor w/ 2111 levels "1 DE MARZO","12 DE DICIEMBRE",..: 553 558 558 2099 2099 1502 1502 346 1611 457 ...
 $ ESQUINA  : Factor w/ 8747 levels "1 DE MARZO Y DUNANT",..: 4372 3657 3657 1630 1630 643 643 6844 7263 5715 ...
 $ geometry :sfc_POINT of length 41121; first list element:  'XY' num  575552 6145117
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA
  ..- attr(*, "names")= chr  "CALLE" "CRUCE" "TIPO" "ANIO" ...

summary(Accidents)

     CALLE          CRUCE         TIPO            ANIO                         NOM_CALLE    
 Min.   :   6   Min.   :   0   FATAL:  543   Min.   :2006   AV ITALIA               : 1417  
 1st Qu.:2052   1st Qu.:2061   GRAVE: 5102   1st Qu.:2007   BV GRAL ARTIGAS         : 1322  
 Median :3639   Median :3780   LEVE :35476   Median :2008   AV GRAL FLORES          : 1143  
 Mean   :3860   Mean   :3867                 Mean   :2008   BV JOSE BATLLE Y ORDOÑEZ:  973  
 3rd Qu.:5880   3rd Qu.:5877                 3rd Qu.:2010   AV GRAL RIVERA          :  967  
 Max.   :9819   Max.   :9819                 Max.   :2010   AV JOSE BELLONI         :  819  
                                                            (Other)                 :34480  
                      NOM_CRUCE                                          ESQUINA         geometry    
 BV JOSE BATLLE Y ORDOÑEZ  :  479   AV ITALIA Y AV CENTENARIO                :  120   POINT  :41121  
 BV GRAL ARTIGAS           :  422   URUGUAYANA Y FRANCISCO GOMEZ             :   87   epsg:NA:    0  
 AV LUIS ALBERTO DE HERRERA:  382   AV ITALIA Y AV BOLIVIA                   :   83                  
 BV APARICIO SARAVIA       :  340   AV GRAL FLORES Y JORGE ISAAC             :   74                  
 AV GRAL FLORES            :  290   AV GRAL FLORES Y BV JOSE BATLLE Y ORDOÑEZ:   67                  
 AV ITALIA                 :  272   AV GRAL RONDEAU Y CNEL FRANCISCO TAJES   :   67                  
 (Other)                   :38936   (Other)                                  :40623                  

By default, plot() creates map for every variable within the file. Use st_geometry to get only the points.

You can add some settings to the maps. Nevetheless, I personally prefer to use tmap.

---

2 Projections (CRS)

Planet earth is not a flat surface and is not perfectly round. Nevertheless, when you make maps you have to put a certain geographic area in two dimensions. This gets even more complicated if you consider mountains and different sea levels.

Projections are the way we transform the shape of the earth (or an area) to only two dimensions. There are many projections and all have some level of distortion. You can see the section 2.4 from the Geocomputation With R book for a clear introduction to the topic. So far, you have to know that every sf object has a Coordinate Referense System (CRS) that can be configured using EPSG codes. There are geographic projections where the notion of distance is not related to a regular understanding of distance (it is not in meter or icnhes). In Projected projections on the other hand, the distance is in metes. Therefore, when you are performing operations sucha as calulations buffers or measuring distances it makes sense to put you CRS as a projected projection.

What is the projection of our Accidents object?

st_crs(Accidents)

Coordinate Reference System: NA

It does not have a defined CRS. You can do that with the st_set_crs() function.

#st_set_crs
#AllCRS <- rgdal::make_EPSG()
#UTM 21 SUR
#32721
#https://epsg.io/32721
#https://www.spatialreference.org/ref/epsg/wgs-84-utm-zone-21s/
Accidents = st_set_crs(Accidents, 32721)
st_crs(Accidents)

Coordinate Reference System:
  EPSG: 32721 
  proj4string: "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

Another usefl function when working with CRS is st_transform than enables to transform the coordinates of an object to other CRS, and switch among projected and geographic references. Be careful of using st_set_crs instead of st_transform when reprojecting and object.

2.1 Tidyverse in action

One of the nicest feautures of sf is that is integrated to the Tidyverse, meaning that you can use the %>% operation and chain your workflow calling dplyr functions such as filter(), slice() or summarize(). An sf object is basically a tibble with an additional column for the geometry.

Accidents %>% names()

[1] "CALLE"     "CRUCE"     "TIPO"      "ANIO"      "NOM_CALLE" "NOM_CRUCE" "ESQUINA"   "geometry" 

Accidents %>% group_by(ANIO) %>% summarise(Total = n())

Simple feature collection with 5 features and 2 fields
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: 558003.5 ymin: 6134508 xmax: 588206.1 ymax: 6157857
epsg (SRID):    32721
proj4string:    +proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs

AccidentsData <- Accidents
st_geometry(AccidentsData) <- NULL
AccidentsData %>% group_by(ANIO) %>% summarise(Total = n())

AccidentsData <- Accidents
AccidentsData %>% group_by(ANIO) %>% summarise(Total = n()) %>% st_set_geometry(NULL)

Resumen <- AccidentsData %>% group_by(ANIO) %>%
  summarise(Total = n()) %>% 
  st_set_geometry(NULL)

ggplot(data = Resumen) + geom_bar(aes(x = ANIO, y = Total), stat = "identity") +
  xlab("Años") + ylab("Accidentes") + theme_classic()

Esquinas <- Accidents %>% group_by(ESQUINA) %>% 
  summarise(Total = n()) %>% 
  st_set_geometry(NULL) %>% 
  arrange(desc(Total))

##Improve this plot!
Esquinas %>% filter(Total >60) %>%  
  ggplot() + 
    geom_bar(aes(x=ESQUINA, y = Total), stat = "identity")

Accidentes2007 <- Accidents %>% filter(ANIO == "2007")
Accidentes2010 <- Accidents %>% filter(ANIO == "2010")
Graves <- Accidents %>% filter(TIPO == "GRAVE")
Graves2008 <- Accidents %>% filter(TIPO == "GRAVE" & ANIO == "2008")
Graves2006o2009 <- Accidents %>% filter(ANIO %in% c("2006","2009"))

---

3 Polygons (States)

With the **st_read()* function is possible to read any kind o vectorial file.

setwd("~/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/Vectoriales_2011")

The working directory was changed to /Users/sabogal/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/Vectoriales_2011 inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

Uruguay <- st_read("ine_depto.shp")

Reading layer `ine_depto' from data source `/Users/sabogal/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/Vectoriales_2011/ine_depto.shp' using driver `ESRI Shapefile'
Simple feature collection with 20 features and 5 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
epsg (SRID):    NA
proj4string:    NA

plot(st_geometry(Uruguay))

Montevideo <- Uruguay[Uruguay$NOMBRE == "MONTEVIDEO",]
plot(st_geometry(Montevideo))

Montevideo <- Uruguay %>% filter(NOMBRE == "MONTEVIDEO")
plot(st_geometry(Montevideo))

st_crs(Montevideo)

Coordinate Reference System: NA

st_crs(Uruguay)

Coordinate Reference System: NA

st_crs(Montevideo) == st_crs(Accidents)

[1] FALSE

Montevideo <- st_set_crs(Montevideo,32721)
st_crs(Montevideo) == st_crs(Accidents)

[1] TRUE

plot(st_geometry(Montevideo))
plot(st_geometry(Accidents), add = TRUE)

---

3.1 The Tidyverse again

Uruguay %>% select(AREA_KM2_) %>% summary()

   AREA_KM2_             geometry 
 Min.   :  248   MULTIPOLYGON:20  
 1st Qu.: 5106   epsg:NA     : 0  
 Median : 9450                    
 Mean   : 8763                    
 3rd Qu.:11714                    
 Max.   :15438                    

Uruguay %>% 
  ggplot(aes(x=AREA_KM2_)) +
    geom_histogram(breaks=seq(0,15000,1000))

BiggesStates <- Uruguay %>% filter(AREA_KM2_ > 13000)
BiggesStates

Simple feature collection with 4 features and 5 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 387005.9 ymin: 6343912 xmax: 858252.1 ymax: 6598411
epsg (SRID):    NA
proj4string:    NA
  AREA_KM2_ PERIMETER DEPTO      NOMBRE CDEPTO_ISO                       geometry
1     13922  712419.7    11    PAYSANDU       UYPA MULTIPOLYGON (((571956.7 64...
2     14163  933919.4    15       SALTO       UYSA MULTIPOLYGON (((418047.2 65...
3     15438  863024.7    18  TACUAREMBO       UYTA MULTIPOLYGON (((565477.8 64...
4     13648  887309.8     4 CERRO LARGO       UYCL MULTIPOLYGON (((694640.5 64...

---

4 Points in Polygons

A simple operation might be to know how many point of an object fall in an area. Let’s use the stations of the railway system in Urugay and see how many stations are in each state. This kind of procedures are called spatial operations.

setwd("~/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/afe_estaciones")

The working directory was changed to /Users/sabogal/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/afe_estaciones inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

Stations <- st_read("afe_estaciones.shp")

Reading layer `afe_estaciones' from data source `/Users/sabogal/Dropbox/Teaching/SpatialAnalysis-MontevideoWorkshop2019/Notebooks/Montevideo_Data/afe_estaciones/afe_estaciones.shp' using driver `ESRI Shapefile'
Simple feature collection with 105 features and 13 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -58.07525 ymin: -34.89368 xmax: -53.38024 ymax: -30.90302
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

st_crs(Stations)

Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Before perfoming any spatial operation is better to ensure that the objects are in the same CRS.

st_crs(Stations) == st_crs(Uruguay)

[1] FALSE

Stations <- st_transform(Stations, 32721)
Uruguay <- st_set_crs(Uruguay, 32721)
st_crs(Stations) == st_crs(Uruguay)

[1] TRUE

plot(st_geometry(Uruguay))
plot(st_geometry(Stations), col = "green",add = T)

StationsInStates <- st_join(Stations, Uruguay)

dim(Stations)

[1] 105  14

dim(Uruguay)

[1] 20  6

dim(StationsInStates)

[1] 105  19

names(Stations)

 [1] "id"         "estado"     "nom_linea"  "tipo"       "nom_estaci" "km"         "cant_desvi"
 [8] "desv_habil" "largo_util" "observacio" "linea_empa" "descriptor" "servicio"   "geometry"  

names(Uruguay)

[1] "AREA_KM2_"  "PERIMETER"  "DEPTO"      "NOMBRE"     "CDEPTO_ISO" "geometry"  

names(StationsInStates)

 [1] "id"         "estado"     "nom_linea"  "tipo"       "nom_estaci" "km"         "cant_desvi"
 [8] "desv_habil" "largo_util" "observacio" "linea_empa" "descriptor" "servicio"   "AREA_KM2_" 
[15] "PERIMETER"  "DEPTO"      "NOMBRE"     "CDEPTO_ISO" "geometry"  

The st_join() funtion performed an operation that retrieved, for every point, the corresponding polygon and attached all data from that polygon. So, to get a polygons object with the amount of points wihtin each state we have to aggregate points by states createing a tibble and join it to the polygons file.

StatesPoints <- StationsInStates %>% group_by(NOMBRE) %>%
  summarise(Stations = n()) %>% 
  st_set_geometry(NULL)
UruguayStations <- Uruguay %>% left_join(StatesPoints)

Joining, by = "NOMBRE"

plot(UruguayStations[,6])

---

5 A better map

So far we have been using the plot() function to create maps. plot() makes the homework but it feels quite limited. I prefer the tmap library and later on the workshop we’ll cover a brief tmap Tutorial. For now, I just want you know that you have the possibility to create atonishing and beautiful maps. But pior to that, there is always a bunch of data processing and wragling that we can do with the help os sf and tidyverse.

tm_shape(UruguayStations) + 
  tm_polygons(col = "Stations",  lwd = 3, border.col = "black") +
  tm_compass(type = "8star", position = c(0.85, 0.04),size=3) +
  tm_scale_bar(position = c(0.06,0))+
  tm_layout(legend.position = c(0.8,0.7),
            legend.title.size = 1,
            legend.frame = T)

---

6 Spatial Operations

6.1 Joining polygons by subsetting

Colonia <- Uruguay %>% filter(NOMBRE == "COLONIA")
Maldonado <- Uruguay %>% filter(NOMBRE == "MALDONADO")

plot(st_geometry(Colonia))

plot(st_geometry(Maldonado))

Union <- rbind(Maldonado, Colonia)
plot(st_geometry(Union))

Union <- rbind(Maldonado, Colonia, Montevideo)
plot(st_geometry(Union))

Montevideo <- Montevideo %>% select(-AREA_KM2_, -NOMBRE)
#Union <- rbind(Maldonado, Colonia, Montevideo) should return an error now. 
Union <- rbind(Maldonado, Colonia)
Union[3,] <- c(rep(0,5),Montevideo[1,4]) #It should be a more elegant way.

invalid factor level, NA generatedinvalid factor level, NA generated

plot(st_geometry(Union))

6.2 Union (dissolve)

#?nc
nc = st_read(system.file("shape/nc.shp", package="sf"))

Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/3.5/Resources/library/sf/shape/nc.shp' using driver `ESRI Shapefile'
Simple feature collection with 100 features and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID):    4267
proj4string:    +proj=longlat +datum=NAD27 +no_defs

plot(st_geometry(nc))

plot(st_union(nc))

nc <- nc %>% mutate(State = "NorthCarolina")
OtherWay <- nc %>%  group_by("NorthCarolina") %>% summarise(Total = n())
plot(st_geometry(OtherWay))

6.3 Buffers

Point1 <- Accidents %>% slice(123)
Point2 <- Accidents %>% slice(321)
#see: https://epsg.io/5383
#https://epsg.io/?q=Uruguay%20kind%3APROJCRS
#Or should we use 54032 ?
Point1 <- st_transform(Point1, 5383)
Point2 <- st_transform(Point2, 5383)
st_is_longlat(Point1)

[1] FALSE

st_is_longlat(Point2)

[1] FALSE

Buffer1 <- st_buffer(Point1, 500)  #The units are un meters
Buffer2 <- st_buffer(Point2, 500)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(Buffer1 %>% st_transform(32721)), add = TRUE)

plot(st_geometry(Buffer2 %>% st_transform(32721)), add = TRUE)

Buffers <- rbind(Buffer1, Buffer2)
Buffers <- st_transform(Buffers, 32721)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(Buffers), add = TRUE)

Points <- Accidents %>% slice(1,2,3,123,321) %>% transform(5383)
BufferPoints <- Points %>% st_buffer(1500) %>% st_transform(32721)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(BufferPoints), col = "gold4",add = TRUE)

plot(st_geometry(Points %>% st_transform(32721)), col = "red", add = TRUE)

6.4 Intersecions

Point1 <- Accidents %>% slice(123)
Point2 <- Accidents %>% slice(321)
Point1 <- st_transform(Point1, 5383)
Point2 <- st_transform(Point2, 5383)
Polygon1 <- st_buffer(Point1,6000) %>% st_transform(32721)
Polygon2 <- st_buffer(Point2,6000) %>% st_transform(32721)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(Polygon1), add = TRUE)

plot(st_geometry(Polygon2), add = TRUE)

MyIntersection <- st_intersection(Polygon1,Polygon2)

attribute variables are assumed to be spatially constant throughout all geometries

#Try st_intersects(Polygon1, Polygon2)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(MyIntersection), col = "darkgreen", border = "darkgreen",add = TRUE)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(Polygon1), add = TRUE)

plot(st_geometry(Polygon2), add = TRUE)
plot(st_geometry(MyIntersection), col = "darkgreen", border = "darkgreen",add = TRUE)

6.5 Difference

MyDifference1 <- st_difference(Polygon1, Polygon2)

attribute variables are assumed to be spatially constant throughout all geometries

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(MyDifference1), col = "darkgreen", border = "darkgreen",add = TRUE)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(Polygon1), add = TRUE)

plot(st_geometry(Polygon2), add = TRUE)
plot(st_geometry(MyDifference1), col = "darkgreen", border = "darkgreen",add = TRUE)

MyDifference2 <- st_difference(Polygon2, Polygon1)

attribute variables are assumed to be spatially constant throughout all geometries

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(MyDifference2), col = "darkgreen", border = "darkgreen",add = TRUE)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(Polygon1), add = TRUE)

plot(st_geometry(Polygon2), add = TRUE)
plot(st_geometry(MyDifference2), col = "darkgreen", border = "darkgreen",add = TRUE)

MyTwoDiffereces <- rbind(MyDifference1, MyDifference2)

plot(st_geometry(Montevideo), col = "deepskyblue1")
plot(st_geometry(Polygon1), add = TRUE)

plot(st_geometry(Polygon2), add = TRUE)
plot(st_geometry(MyTwoDiffereces), col = "darkgreen", border = "darkgreen",add = TRUE)

6.6 Areas

st_area(Polygon1)

112457788 [m^2]

st_area(Polygon2)

112449284 [m^2]

Uruguay <- Uruguay %>% mutate(Areas = st_area(Uruguay)) %>% mutate(Areas = as.numeric(Areas))

6.7 Ceontroids

MyCentroid <- st_centroid(Polygon1)

st_centroid assumes attributes are constant over geometries of x

MyCentroid

Simple feature collection with 1 feature and 7 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 581305.4 ymin: 6139145 xmax: 581305.4 ymax: 6139145
epsg (SRID):    32721
proj4string:    +proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs
  CALLE CRUCE TIPO ANIO NOM_CALLE         NOM_CRUCE                       ESQUINA
1  3564  3525 LEVE 2010 AV ITALIA HIPOLITO YRIGOYEN AV ITALIA Y HIPOLITO YRIGOYEN
                  geometry
1 POINT (581305.4 6139145)

plot(st_geometry(Polygon1), col = "black")
plot(st_geometry(MyCentroid), col = "red", pch = 3, cex = 4, add = TRUE)

UruguayCentroids <- Uruguay %>% st_centroid()

st_centroid assumes attributes are constant over geometries of x

plot(st_geometry(Uruguay), col = "deepskyblue1")
plot(st_geometry(UruguayCentroids), pch = 3, col = "red", add = TRUE)

6.8 Other Spatial Operations

If you want to see the full capabilities provided by sf you can check this cheatsheet or the package vignettes

---

7 Data Sources

All data used in this workflow was downloaded from Open Source pltaforms

• Data in /montevideorural_shp and /montevideourbano_shp were retrieved from the Uruguay Cátalogo dde datos abiertos.
• Data in /Vectoriales_2011 were retrieved from Instituto nacional de estadística
• Data in /Accidentes2006-2010 were taken from Catálogo de datos abiertos
• Data in /CyclingMontevideoData were taken from Catálogo de datos abiertos
• Data in/afe_estaciones were taken from the Geoportal

Another interesting soure in Uruguay is:

• La Encuesta continua de hogares