Tarea 2. Simple Features for R

En esta tarea vamos a representar las pistas de tenis de la Comunidad Valenciana por províncias.

Cargamos las librerías

library(tidyverse)
library(sf)
library(osmdata)
library(ggmap)

Cargamos los Datos

Vamos a cargar datos de Cartografía digitalizada de secciones censales en formato shp,descargados del Instituto Nacional de Estadística. Posteriormente, obtenemos aquellos que se corresponden con la Comunidad Valenciana.

data <- st_read("Seccionado_2021/SECC_CE_20210101.shp")

## Reading layer `SECC_CE_20210101' from data source 
##   `C:\Users\sanch\OneDrive\Escritorio\BIA\2021-2022\Datos_Espaciales_y_Espacio_Temporales\practica3\Seccionado_2021\SECC_CE_20210101.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 36333 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -1004502 ymin: 3132130 xmax: 1126932 ymax: 4859240
## Projected CRS: ETRS89 / UTM zone 30N

CV_censal<-subset(data, NCA=="Comunitat Valenciana")

Representamos gráficamente:

plot(st_geometry(CV_censal))

Se muestra la Comunidad Valenciana dividida por censos. Por lo que para obtener el mapa de la Comunidad Valenciana “limpio” agrupamos los censos por la comunidad autonoma a la que pertenecen:

CV<- CV_censal%>%
  group_by(CCA)%>%
  summarise()
plot(st_geometry(CV))

En nuestro caso, vamos a representar las pistas de tenis en la Comunidad Valenciana pero diferenciando cada provincia.Para obtener el mapa de la Comunidad Valenciana dividido por provincias agrupamos los censos por la provincia a la que pertenecen:

CV_Prov<- CV_censal%>%
  group_by(CPRO) %>%
  summarise()
plot(st_geometry(CV_Prov))

Como bien hemos comentado, queremos representar geográficamente las pistas de tenis de la Comunidad Valenciana por provincias. Obtenemos los datos de la web de Wiki Open Street Map (https://wiki.openstreetmap.org/wiki/Map_features) y usaremos la etiqueta de “tennis” de la característica “sport”. Creamos la capa de puntos:

assign("has_internet_via_proxy", TRUE, environment(curl::has_internet))
m <- c(-1.5,37.5,0.5,41) # Oeste/Sur/Este/Norte
q1 <- m %>%
  opq(timeout = 2500) %>%
  add_osm_feature("sport", "tennis")

tennis <- osmdata_sf(q1)

ggplot(tennis$osm_points) +
  geom_sf(colour = "#339900",
          fill = "#66FF00",
          alpha = .5,
          size = 1,
          shape = 21) +
  theme_void()

Ahora tenemos el mapa por un lado y los puntos por otro. El siguiente paso consiste en unir las dos capas. Unimos el mapa con los puntos y representamos:

ggplot(CV_Prov)+
  geom_sf()+ theme_bw()+
  geom_sf(fill = "white", color = "black")+
  xlab("Longitud") + ylab("Latitud") +
  geom_sf(data = tennis$osm_points,
          inherit.aes = FALSE,
          colour = "#339900",
          fill = "#66FF00",
          alpha = .5,
          size = 1,
          shape = 21) +
  labs(x = "", y = "")+
  scale_color_brewer(type = 'div', palette = 'Set1', direction = 1)+
  theme_minimal()

Vemos que hay puntos que no pertenecen a la Comunidad Valenciana, procedemos a eliminarlos. Tenemos el dataset CV_Prov donde tenemos polígonos con los límites de cada provincia de la Comunidad Valenciana. Por otro lado datos contiene los puntos de las pistas de tenis. Podemos usar la función st_join que devuelve un nuevo objeto sf con los atributos que corresponden a la intersección espacial entre los puntos y polígonos.

datos<-tennis$osm_points
# datosTennisCV <- st_join(datos, CV_Prov)

Nos devuelve un error ya que los datos están en CRS distintos.

st_crs(datos)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

st_crs(CV_Prov)

## Coordinate Reference System:
##   User input: ETRS89 / UTM zone 30N 
##   wkt:
## PROJCRS["ETRS89 / UTM zone 30N",
##     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["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe between 6Â°W and 0Â°W: Faroe Islands offshore; Ireland - offshore; Jan Mayen - offshore; Norway including Svalbard - offshore; Spain - onshore and offshore."],
##         BBOX[35.26,-6,80.53,0]],
##     ID["EPSG",25830]]

Hacemos la transformación con la función st_transform y después hacemos el join.

datos2 <- datos %>% st_transform(st_crs(CV_Prov))
datosTennisCV <- st_join(x = datos2, y = CV_Prov)
colnames(datosTennisCV)

##  [1] "osm_id"           "name"             "access"           "addr.housenumber"
##  [5] "addr.street"      "amenity"          "barrier"          "building"        
##  [9] "created_by"       "entrance"         "leisure"          "nohousenumber"   
## [13] "source"           "sport"            "CPRO"             "geometry"

Creamos una nueva variable que diga TRUE si hay algún dato en esa variable y FALSE cuando no.

datosTennisCV <- datosTennisCV %>%
  mutate(CABA=ifelse(is.na(CPRO),TRUE,FALSE))

table(datosTennisCV$CABA)

## 
## FALSE  TRUE 
##  9952  2036

De las 11705 observaciones, 1879 no se corresponden con nuestro polígono del mapa de la Comunidad Valenciana. Podemos representar graficamente las pistas de tenis que están en la Comunidad Valenciana de color rojo y de color azul las que se encuentren fuera de esta comunidad autónoma.

ggplot(CV_Prov)+
  geom_sf()+ theme_bw()+
  geom_sf(fill = "white", color = "black")+
  xlab("Longitud") + ylab("Latitud") +
  ggtitle("Mapa de la Comunidad Valenciana")+
  geom_sf(data = datosTennisCV,
          mapping = aes(color=CABA, fill=CABA),
          inherit.aes = FALSE,
          alpha = .8,
          size = 1,
          shape = 21,
          show.legend = FALSE) +
  labs(x = "", y = "")+
  scale_color_brewer(type = 'div', palette = 'Set1', direction = 1)+
  theme_minimal()

Nos interesan los puntos que figuran en color rojo, eliminamos por tanto los de color azul Para ello hacemos uso de la función filter.

datos_definitivos<- datosTennisCV %>%
                      filter(CABA == FALSE)


ggplot(CV_Prov)+
  geom_sf()+ theme_bw()+
  geom_sf(fill = "white", color = "black")+
  xlab("Longitud") + ylab("Latitud") +
  ggtitle("Mapa de la Comunidad Valenciana")+
  geom_sf(data = datos_definitivos,
          mapping = aes(color="#99FF33", fill="#99FF33"),
          inherit.aes = FALSE,
          alpha = .8,
          size = 1,
          shape = 21,
          show.legend = FALSE) +
  labs(x = "", y = "")+
  scale_color_brewer(type = 'div', palette = 'Set1', direction = 1)+
  theme_minimal()

Añadimos contexto geográfico

Para dar un mayor contexto geográfico, añadimos parte del territorio de las comunidades autónomas colindantes a la Comunidad Valenciana. Creamos un objeto que agrupe el mapa censal de españa en comunidades autónomas y creamos un objeto que contenga aquellas comunidades autónomas que tocan las provincias de la Comunidad Valenciana. Finalmente, representamos gráficamente el resultado:

CAMAP<- data%>%
  group_by(CCA)%>%
  summarise()
toca_CV<-CAMAP[CV_Prov, op = st_touches]
plot(st_geometry(toca_CV),col= "#999999")

Finalemnte, lo añadímos a nuestro mapa y mostramos el mapa final.

limsCV <- st_buffer(CV_Prov,dist = 0.3) %>% st_bbox()

ggplot()+
  geom_sf(data = toca_CV, colour="#000000", fill= "#999999") +
  geom_sf(data = CV_Prov, colour="#000000", fill = "#CCCCCC") +
  geom_sf(data = datos_definitivos,
          inherit.aes = FALSE,
          colour= "#339900", 
          fill="#66FF00",  
          alpha = .8,
          size = 1.3,
          shape = 21,
          show.legend = FALSE) +
  scale_alpha_continuous(range = c(.6, 1), guide=FALSE) +
  scale_size_identity(guide = FALSE) +
  coord_sf(xlim = c(limsCV["xmin"], limsCV["xmax"]), 
           ylim = c(limsCV["ymin"], limsCV["ymax"])) +
  xlab("Longitud") + ylab("Latitud") +
  labs(title = "Pistas de tenis en la Comunidad Valenciana \n por provincias.",
       subtitle = "Realizado por: Marina Sanchis") +
    theme(panel.grid = element_line(colour = 'transparent'),
        line = element_blank(), 
        rect = element_blank())+
  theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", 
                                        size = 0.5),
        panel.background = element_rect(fill = "aliceblue"))+
  scale_color_brewer(type = 'div', palette = 'Set1', direction = 1)

## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.
## It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

Tarea 2. Simple Features for R

Marina Sanchis Romá

27/2/2022

Cargamos las librerías

Cargamos los Datos

Añadimos contexto geográfico