I'm new to programming in R and I want to make an interactive map from two files, one is a .shp that you can download from here: https://www.ine.es/ss/Satellite?L=es_ES&c=Page&cid=1259952026632&p=1259952026632&pagename=ProductosYServicios%2FPYSLayout (just select 2021 year and go and its download), in which there are many polygons. And then I have a csv with store characterization data (it contains 2 LON and LAT fields).
To start doing all this I would like to filter the .shp file for each different value in the NCA field (Ex: 1 map for Basque Country, another for Madrid, another for Barcelona ...).
All this without losing the geometric properties since if I lose them then I can't represent them graphically (or maybe I can and I don't know, if so, let me know and I will be very grateful).
He probado con el siguiente codigo:
# Load the libraries
pacman::p_load(leaflet, leaflet.extras, mapview, rworldxtra, rgdal,raster, sf, tidyverse, readr, ggthemes)
# Load the .shp file in spdf format.
myspdf = readOGR(getwd(), layer = "SECC_CE_20210101")
#Filter
PV = myspdf %>% filter(NCA == "País Vasco") # Dont work
PV2 = myspdf[myspdf$NCA == "País Vasco"] # Dont work
When I load the shp file and save it in the variable myspdf, I can visualize something like this: https://ibb.co/mywDd6p
in which if I do myspdf@data I access the data (where is the NCA field where I want to filter)
So when I try to filter like this:
PV = myspdf %>% filter(NCA == "País Vasco") # Dont work
PV2 = myspdf[myspdf$NCA == "País Vasco"] # Dont work
It returns this to me this: https://ibb.co/VDYdByq, with the rows completely empty, and what I would like to obtain is the same format but with about 1700 rows x 18 columns and with the geometric properties as well.
Another question I have is that when I read the .shp file as sf, one more column is added with the geometry and inside are the coordinates stored in lists, like that: https://ibb.co/M1Fn8K5, I can easily filter it but I don't know how to represent it graphically (leaflet or mapview...) so that You can see the polygons of NCA = 'Basque Country', could someone give me an example with this?
Ok! I guess I will do the all workflow!
library(sf)
library(tmap)
library(mapview)
# lets get some shops
shop <- data.frame(X = c(-4.758628, -4.758244, -4.756829, -4.759394, -4.753698,
-4.735330, -4.864548, -4.863816, -4.784694, -4.738924),
Y = c(43.42144, 43.42244, 43.42063, 43.42170, 43.41899,
43.41181, 43.42327, 43.42370, 43.42422, 43.40150),
name = LETTERS[1:10])
# Here I save them
write.csv(shop, "shop.csv")
# because I want to show you how to import
shop <- read.csv("shop.csv")
# and convert to en sf object
shop_sf <- sf::st_as_sf(shop, coords = c("X", "Y"))
# and add a CRS
shop_sf <- sf::st_set_crs(shop_sf, 4326)
# now I have downloaded data from your link
# I import it in R
spain_seccionado <- sf::st_read("España_Seccionado2021/SECC_CE_20210101.shp")
# Rq CRS is ETRS89 / UTM 30, will need to transform that
# here I am just exploring a bit the data set
names(spain_seccionado)
unique(spain_seccionado$NCA)
# I just keep Asturias, You have plenty of different way of doing that
# this is what you tried to do here: PV = myspdf %>% filter(NCA == "País Vasco")
# but on an sp object not an sf one
Asturias <- spain_seccionado[spain_seccionado$NCA == "Principado de Asturias",]
asturias_4326 <- sf::st_transform(Asturias, 4326)
# Now both data set are in the same CRS
# a quick plot just to see if everything is correct
plot(asturias_4326$geometry)
plot(shop_sf, col = "red", add = TRUE, pch = 5)
# An interactive map quick and dirty you will need to improve it !
tmap_mode("view")
llanes_shop <- tmap::tm_shape(asturias_4326) +
tmap::tm_borders() +
tmap::tm_shape(shop_sf) +
tmap::tm_symbols(shape = 24) +
tmap::tm_layout()
llanes_shop