Accidentesyvialidades
Claudia Guerrero S
2025-10-11
Accidentes de Tránsito terrestrestre en zonas urbanas y suburbanas (ATUS) https://www.inegi.org.mx/programas/accidentes/#microdatos
#BASE MUNICIPAL_ACCIDENTES DE TRANSITO GEORREFERENCIADOS_2023
ruta<-"../../data/atus_2023_shp/conjunto_de_datos/"
infile <- "BASE MUNICIPAL_ACCIDENTES DE TRANSITO GEORREFERENCIADOS_2023.shp"
nomarchi<-as.filename(paste0(ruta,infile))
nomfil<-"BASE MUNICIPAL_ACCIDENTES DE TRANSITO GEORREFERENCIADOS_2023"
datacom<- read_sf(nomarchi,options = "ENCODING=latin1")
datacom<-st_transform(datacom, CRS("+init=epsg:4326 +proj=longlat
+ +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
datacom_sub <-subset(datacom,(datacom$EDO %in% c(5) & datacom$MPIO %in% c(35)))
datacom_sub <-datacom_sub |>
filter(MOTOCICLET >0)Voronoi
#Voronoi no puntos duplicados
a<-as(datacom_sub, "Spatial")
class(a)[1] “SpatialPointsDataFrame” attr(,“package”) [1] “sp”
min_dist<-0.1
datacomunique<-sp::remove.duplicates(a, zero = min_dist)
class(datacomunique)[1] “SpatialPointsDataFrame” attr(,“package”) [1] “sp”
datacomuniquesf<- st_as_sf(datacomunique)
class(datacomuniquesf)[1] “sf” “data.frame”
#### para voronoi solo un punto
df.2 <- datacomuniquesf[, c("LONGITUD","LATITUD","geometry")]
colnames(df.2)[1]<-"x"
colnames(df.2)[2]<-"y"
#remotes::install_github("garretrc/ggvoronoi", dependencies = TRUE, build_opts = c("--no-resave-data"))
ggplot(df.2 , aes(x, y)) +
stat_voronoi(geom = "path")+
geom_point(aes(x=x,y=y))
###BASE mapa geoestadística nacional
#nomfile<-'C:/users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/Marcogeoestadisticonacional/mg_2023_integrado/conjunto_de_datos/00ent.shp'
#nomfil<-"00ent.shp"
#ent<-read_sf(nomfile, options = "ENCODING=latin1")
#ent<-st_transform(ent, CRS("+init=epsg:4326 +proj=longlat
#+ +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
#st_write(ent, 'ent.gpkg', append=FALSE)
ruta<-"C:/users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/Marcogeoestadisticonacional/mg_2023_integrado/conjunto_de_datos/"
infile <- "00mun.shp"
nomarchi<-as.filename(paste0(ruta,infile))
nomfil<-"00mun"
a<-"SELECT * FROM \"00mun\" WHERE CVE_ENT = '05'"
#a<-"SELECT * FROM \"00mun\" WHERE CVEGEO = '05017' OR CVEGEO = '05033' OR CVEGEO = '05035' OR CVEGEO = '10007' OR CVEGEO = '10012'"
datageo<- read_sf(nomarchi,options = "ENCODING=UTF-8",query=a)
#datageo<- read_sf(nomarchi,options = "ENCODING=UTF-8")
datageo<-st_transform(datageo, CRS("+init=epsg:4326 +proj=longlat
+ +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
#st_geometry(datageo)
datageo_sub <-subset(datageo,datageo$CVE_MUN %in% c("035"))
#datageo_sub$_ogr_geometry_
attr(datageo_sub, "sf_column")[1] “ogr_geometry”
datageo_sub_geom <- st_geometry(datageo_sub)
polygon = st_cast(datageo_sub_geom, "POLYGON")
#plot(polygon)
tor<-polygon[2]
#plot(tor)
#bb<-st_as_text(st_geometry(datageo_sub))
#bb<-st_as_text(st_geometry(polygon))
#bb<-st_as_text(st_geometry(polygon))
#### mapa de red carreteras inegi
viasredvial<- read_sf('C:/users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/Marcogeoestadisticonacional/Carreteras/conjunto_de_datos/red_vial.shp')
viasredvial<-st_transform(viasredvial, CRS("+init=epsg:4326 +proj=longlat
+ +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
viasredvial_sub <-viasredvial[tor,]
pal <- "monet"
colors<-met.brewer("Monet")
#Monet = list(c("#4e6d58", "#749e89", "#abccbe", "#e3cacf", "#c399a2", "#9f6e71", "#41507b", "#7d87b2", "#c2cae3"), c(2, 5, 8, 3, 4, 9, 1, 6, 7), colorblind=FALSE),
#swatchplot(colors)
texture <- grDevices::colorRampPalette(colors, bias = 4)(256)
#swatchplot(texture)
theme_map <- function(...) {
theme_minimal() +
theme(
legend.position="bottom",
legend.direction="horizontal",
#legend.key.size = unit(2, "cm"),
#legend.text=element_text(size=15),
text = element_text(family = "Sans Serif", color = "#c2cae3"),
strip.text = element_text(face = "bold", size = rel(1.5)),
strip.background = element_rect(fill = colors[1], colour = "black"),
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
#panel.grid.minor = element_line(color = colors[8], size = 0.2),
panel.grid.major = element_line(color = colors[2], linewidth = 1),
panel.grid.minor = element_blank(),
#panel.grid.major = element_blank(),
plot.background = element_rect(fill = colors[11], color = "transparent"),
#plot.background = element_blank(),
panel.background = element_rect(fill = colors[11], color = "transparent"),
#panel.background = element_blank(),
legend.background = element_rect(fill = "transparent", color = "transparent"),
panel.border = element_blank(),
...
)
}
`%notin%` <- Negate(`%in%`)
#unique(viasredvial_sub$CARRILES)
#unique(recorteviasredvial_sub$CARRILES)
#mapamx<-ggplot(polygon[2])+
mapamx<-ggplot(tor)+
geom_sf(fill = colors[9],color=colors[7],size=2,alpha=.1) +
geom_sf(data = datacom_sub,color=colors[7],size=2,alpha=1)+
theme_map() +
coord_sf()
mapamx
#rtistry280
mapamx<-ggplot(df.2,aes(x,y))+
stat_voronoi(geom = "path",
color = colors[2], # Color of the lines
lwd = 0.5, # Width of the lines
linetype = 1) + # Type of the lines
geom_sf(data = df.2,color=colors[7],size=1,alpha=1)+
geom_sf(data = tor, inherit.aes = FALSE,color=colors[6],alpha=.10)+
theme_map() +
coord_sf()
mapamx
polygon2sf<- st_as_sf(polygon[2])
#class(polygon2sf)
dft <- polygon2sf %>%
# 2 Extract coordinates
st_coordinates() %>%
# 3 to table /tibble
as.data.frame()
dft.2 <- dft[, c("X","Y")]
colnames(dft.2)[1]<-"x"
colnames(dft.2)[2]<-"y"
dft.2$group<-1
#swatchplot(colors)
mapamx<-ggplot(df.2,aes(x,y))+
geom_voronoi(outline = dft.2,
color = colors[7], # Color of the lines
lwd = 0.8, # Width of the lines
linetype = 1,# Type of the lines
alpha=0.01) +
geom_sf(data = df.2,color=colors[5],size=.1,alpha=1)+
geom_sf(data = tor, inherit.aes = FALSE,color=colors[6],alpha=.01,lwd = 1)+
theme_map() +
coord_sf()
mapamx
if (!dir.exists(glue("images/{map}"))) {
dir.create(glue("images/{map}"))
}
# Set up outfile where graphic will be saved.
# Note that I am not tracking the `images` directory, and this
# is because these files are big enough to make tracking them on
# GitHub difficult.
outfile <- str_to_lower(glue("images/{map}/{map}_{pal}M2.png"))
# Now that everything is assigned, save these objects so we
# can use then in our markup script
width = NA
height = 800
asp <- get_aspect_ratio(tor)
if (is.na(width) && !is.na(height)) {
width <- height * asp
} else if (is.na(height) && !is.na(width)) {
height <- width / asp
}
#"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/portraits/torreon/images/torreon"
#https://stackoverflow.com/questions/57868862/crop-unused-margin-in-a-ggplot-plot
#ggsave(outfile, width = 800, height = 800, units = "mm", dpi = "retina",bg="white",limitsize = FALSE)
ggsave(outfile, width = width, height = height, units = "mm", dpi = "retina",bg="white",limitsize = FALSE)Hexágonos
sf::sf_use_s2(FALSE)
hexgrid <- st_make_grid(tor,
cellsize = .005, ## unit: metres; change as required
what = 'polygons',
square = FALSE ## !
) |>
st_as_sf()
#plot(hexgrid)
hexgrid_tor <- hexgrid[c(unlist(st_contains(tor, hexgrid)),
unlist(st_overlaps(tor, hexgrid))) ,]
#tor|> plot()
#hexgrid_tor |> plot(add = TRUE)
grid_spacing <- 0.01
poly <- st_make_grid(tor, square = F, what="polygons",
#cellsize = c(grid_spacing, grid_spacing*0.5)) %>%
cellsize = grid_spacing*0.25) %>%
st_intersection(tor) %>%
st_sf()
## plot
#ggplot()+
# geom_sf(data=poly, color="red")+
# theme_void()
#BASE MUNICIPAL_ACCIDENTES DE TRANSITO GEORREFERENCIADOS_2023
## plot
ggplot()+
geom_sf(data=poly,fill = alpha(colors[5], 0.2), color= alpha(colors[6], 1))+
geom_sf(data = datacom_sub,color=colors[8],size=.5,alpha=1)+
theme_void()
#https://urbandatapalette.com/post/2021-08-tessellation-sf/
# To sf and add grid ID
poly_sf = st_sf(poly) %>%
# add grid ID
mutate(grid_id = 1:length(lengths(poly)))
# count number of points in each grid
# https://gis.stackexchange.com/questions/323698/counting-points-in-polygons-with-sf-package-of-r
poly_sf$n_colli = lengths(st_intersects(poly_sf, datacom_sub))
poly_sf_count = filter(poly_sf, n_colli > 0)
#poly_sf_count = filter(poly_sf, n_colli > -1)
#ggplot()+
#geom_sf(data = poly_sf_count,color=colors[1],size=.5,alpha=1)+
# theme_void()
library(cptcity)
pals <- find_cpt("pop")
#http://soliton.vm.bytemark.co.uk/pub/cpt-city/notes/formats.html
#http://soliton.vm.bytemark.co.uk/pub/cpt-city/index.html
pals <- find_cpt("ncl")
pals <- find_cpt("seq")
breaks <- classInt::classIntervals(
poly_sf_count$n_colli,
n = 15,
# style = "jenks"
style = "fisher"
)$brks
#scale_fill_gradientn(colours = cpt(n = 10, "cb_seq_PuBuGn_09"),breaks = round(breaks, 0))+
ggplot()+
#geom_sf(data=poly,fill = alpha(colors[5], 0.2), color= alpha(colors[6], 1))+
geom_sf(data = tor,fill="white",color="black",size=2,alpha=1)+
geom_sf(data = poly_sf_count,aes(fill=n_colli),size=.5,alpha=1)+
scale_fill_gradientn(colours = cpt(n = 15, "jjg_cbcont_seq_cbcPurples"),breaks = round(breaks, 0))+
geom_sf(data = datacom_sub,color="orange",size=.1,alpha=1)+
theme_map() +
coord_sf()
outfile <- str_to_lower(glue("images/{map}/{map}_{pal}M3.png"))
# Now that everything is assigned, save these objects so we
# can use then in our markup script
width = NA
height = 800
asp <- get_aspect_ratio(tor)
if (is.na(width) && !is.na(height)) {
width <- height * asp
} else if (is.na(height) && !is.na(width)) {
height <- width / asp
}
#"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/portraits/torreon/images/torreon"
#https://stackoverflow.com/questions/57868862/crop-unused-margin-in-a-ggplot-plot
#ggsave(outfile, width = 800, height = 800, units = "mm", dpi = "retina",bg="white",limitsize = FALSE)
ggsave(outfile, width = width, height = height, units = "mm", dpi = "retina",bg="white",limitsize = FALSE)Heat
pal <- "Hokusai1"
colors3<-met.brewer("Egypt",n=4)
#swatchplot(colors3)
colors3[4][1] “#fab255”
Areas y zonas elegidas uso de osm
#Creación del objeto espacial seleccionado y su proyección
### recorte de zona
rutacam<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/"
limitesbox<-read.csv(paste0(rutacam,"limites.txt"),header=TRUE, sep=",",encoding="latin")
names(limitesbox)<-c("long","lat")
#Creación del objeto espacial seleccionado y su proyección
sps <- st_as_sf(limitesbox, coords = c("long","lat"),crs = 4326)%>%
summarise(geometry = st_combine(geometry)) %>%
st_cast("POLYGON")
sps<-st_transform(sps, CRS("+init=epsg:4326 +proj=longlat
+ +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
bb <- sps %>%
pull(geometry) %>%
st_as_text()
lim<-st_bbox(sps)
#plot(qosm_sf$geometry)
xlim <- c(lim$xmin,lim$xmax)
ylim <- c(lim$ymin,lim$ymax)
nomarchi<-as.filename("C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/gis_osm_roads_free_1.gpkg")
nomfil<-"gis_osm_roads_free_1"
a<-glue("SELECT * FROM ", nomfil, " WHERE st_intersects(geom, st_polygonfromtext(")
c<-paste0(a,"'",bb,"' ))")
#c<-paste0(a,"'",bb1,"' )) OR st_intersects(geom, st_polygonfromtext(","'",bb2,"' )) ")
#viasosmlayer_sub<-read_sf(nomarchi, query = c)
viasosmlayer_sub<-read_sf(nomarchi)
viasosmlayer_sub<-st_transform(viasosmlayer_sub, CRS("+init=epsg:4326 +proj=longlat
+ +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
viasosmlayer_sub <-viasosmlayer_sub[tor,]
#Ya recortada en la lectura... si no por subconjunto
#recorteviasosmlayer_sub<-viasosmlayer_sub[sps, ]
recorteviasosmlayer_sub<-viasosmlayer_sub
viasredvial_sub<-st_transform(viasredvial_sub, CRS("+init=epsg:4326 +proj=longlat
+ +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
#Ya recortado desde la consulta
#recorteviasredvial_sub<-viasredvial_sub[sps, ]
recorteviasredvial_sub<-viasredvial_sub
#### mapa de red carreteras inegi
#viasredvial_sub<- read_sf('D:/Documents/Claudia/Midropbox/Investigacion y escritos/Marcogeoestadisticonacional/Carreteras/conjunto_de_datos/red_vial.shp')
#BASE MUNICIPAL_ACCIDENTES DE TRANSITO GEORREFERENCIADOS_2023
#este archivo no se leyó con consulta, se recortó después
#recortedatalum <-datacom_sub[sps, ]
recortedatalum <-datacom_sub
#"Clave 0.- Certificado cero
#Clave 1.- Colisión con vehículo automotor
#Clave 2.- Colisión con peatón (atropellamiento)
#Clave 3.- Colisión con animal
#Clave 4.- Colisión con objeto fijo
#Clave 5.- Volcadura
#Clave 6.- Caída de pasajero
#Clave 7.- Salida del camino
#Clave 8.- Incendio
#Clave 9.- Colisión con ferrocarril
#Clave 10.- Colisión con motocicleta
#Clave 11.- Colisión con ciclista
#Clave 12.- Otro"
#unique(datacom_sub$TIPACCID)
recortedatalum <-datacom_sub |>
filter(MOTOCICLET >0)
####recodificación diccionario de datos
dias<-c("Lunes","Martes","Miércoles","Jueves","Viernes","Sábado","Domingo","No especificado")
recortedatalum$DIASEMANA<-dias[as.integer(recortedatalum$DIASEMANA)]
tipac<-c("clave cero","Colisión automotor","atropellamiento","Colisión con animal","Colisión con objeto fijo","Volcadura","Caída de pasajero","Salida del camino","Incendio","Colisión con ferrocarril","Colisión con motocicleta","Colisión con ciclista","Otro")
recortedatalum$TIPACCID<-tipac[as.integer(recortedatalum$TIPACCID)]
sexop<-c("Se fugó","Hombre","Mujer")
recortedatalum$SEXO<-sexop[as.integer(recortedatalum$SEXO)]
grave<-c("Fatal","No fatal","Sólo daños")
recortedatalum$CLASE<-grave[as.integer(recortedatalum$CLASE)]
unique(recortedatalum$CLASE)[1] “No fatal” “Sólo daños” “Fatal”
# set up color palette
c1 <- PrettyCols::prettycols("Greens")
c2 <- PrettyCols::prettycols("Tangerines")
# set up color palette
#https://github.com/BlakeRMills/MetBrewer
pal <- "metmonetgreenbako"
colors<-met.brewer("Monet")
#Monet = list(c("#4e6d58", "#749e89", "#abccbe", "#e3cacf", "#c399a2", "#9f6e71", "#41507b", "#7d87b2", "#c2cae3"), c(2, 5, 8, 3, 4, 9, 1, 6, 7), colorblind=FALSE),
#swatchplot(colors)
texture <- grDevices::colorRampPalette(colors, bias = 4)(256)
#swatchplot(texture)
colors2<-met.brewer("Troy",n=13)
#swatchplot(colors2)
texture <- grDevices::colorRampPalette(colors2, bias = 4)(256)
#swatchplot(texture)
colors3 <- scico::scico(n = 300, palette = "bamako")
#swatchplot(colors3)
texture <- grDevices::colorRampPalette(colors3, bias = 4)(256)
#swatchplot(texture)
#Polígonos de densidad
#CAPAS de polígonos
#http://www.mosaic-web.org/ggformula/articles/pkgdown/ggformula-long.html
#https://gis.stackexchange.com/questions/444689/create-density-polygons-from-spatial-points-in-r
#Calculate the centroid of each county of interest
# https://stackoverflow.com/questions/49343958/do-the-values-returned-by-rgeosgcentroid-and-sfst-centroid-differ
#NatAM_counties <- counties_shp |>
# filter(name_origin == "Native American origin")
#centers se usa como dato de la densidad en el geom_den_sity
#centers <- data.frame(
# st_coordinates(st_centroid(NatAM_counties$geometry)),
# id=NatAM_counties$NAME)
p <- ggplot(recortedatalum, aes(x = purrr::map_dbl(geometry, ~.[1]),
y = purrr::map_dbl(geometry, ~.[2])))+
geom_density2d_filled()+
geom_point() +
geom_density2d()
p
#https://stackoverflow.com/questions/30626400/difference-in-2d-kde-produced-using-kde2d-r-and-ksdensity2d-matlab
#kde2d(x, y, h, n = 25, lims = c(range(x), range(y)))
#f2 <- kde2d(duration, waiting, n = 50, lims = c(0.5, 6, 40, 100),
# h = c(width.SJ(duration), width.SJ(waiting)) )
######Estructura del objeto
#gb <- ggplot_build(p)
#gb$data
#gb$layout
#gb$plot
#gb$data[1]
#gb$data[2]
#gb$data[3]
#unlist(gb$data[[1]])
#gb$data[c(1,2)] #dos elementos de la lista (tiene 3)
#gb$data[[3]][[1]] #el primer elemento de la lista 3 de la lista data
data2d <- layer_data(p, i = 1)
# Id of polygon
data2d$pol <- paste0(data2d$group, "_", data2d$subgroup)
ids <- unique(data2d$pol)
# Split and create polygons based on the id
pols <- lapply(ids, function(x){
topol <- data2d[data2d$pol == x, ]
closepol <- rbind(topol, topol[1, ])
pol <- st_polygon(list(as.matrix(closepol[,c("x", "y")])))
# Add features
df <- unique(topol[, grepl("level", names(topol))])
tofeatures <- st_as_sf(df, geometry=st_sfc(pol))
return(tofeatures)
})
final_pols <- do.call(rbind, pols)
# And force a crs, since we lost that on the process
st_crs(final_pols) <- st_crs(viasosmlayer_sub)
#final_pols$level
nc_geom <- st_geometry(final_pols)
#plot(nc_geom)
label = levels(final_pols$level)
#geometria primer nivel
#filterpols<-filter(final_pols,level %in% c("(0, 50]","(50, 100]","(100, 150]","(150, 200]","(200, 250]","(250, 300]","(300, 350]"))
filterpols<-filter(final_pols,level %in% label)
#Selección de polígonos
filterpols2 <- filterpols[-1, ]
#shp.subset <- shp[i %in% c(1,2,3),]
#shp.sub <- shp[shp$AREA < 10,]
#shp.sub <- subset(shp, AREA < 10)
`%notin%` <- Negate(`%in%`)
#unique(recorteviasredvial_sub$CARRILES)
mapamx<-ggplot()+
geom_sf(data=filterpols2,aes(fill=level),show.legend=TRUE)+
scale_fill_manual(values=colors2) +
ggnewscale::new_scale_color() +
geom_sf(data=viasosmlayer_sub |>
filter(fclass %in% c("primary","primary_link")),
inherit.aes = FALSE,
color = colors[5],
size = 0.8,
alpha=.6)+
geom_sf(data=viasosmlayer_sub |>
filter(fclass %in% c("scondary","tertiary","residential")),
inherit.aes = FALSE,
color = colors[3],
size = 0.6,
alpha=.6)+
geom_sf(data=viasosmlayer_sub |>
filter(fclass %notin% c("primary","primary_link","scondary","tertiary","residential")),
inherit.aes = FALSE,
color = colors[8],
size = 0.5,
alpha=.6)+
geom_sf(data=viasredvial_sub |>
filter(CARRILES %in% c("6","5","4")),
inherit.aes = FALSE,
color = colors[1],
size = 1,
alpha=.6)+
geom_sf(data=viasredvial_sub |>
filter(CARRILES %in% c("3","2","1")),
inherit.aes = FALSE,
color = colors[1],
size = 1,
alpha=.6)+
geom_sf(data=viasredvial_sub |>
filter(CARRILES %notin% c("6","5,","4","3","2","1")),
inherit.aes = FALSE,
color = colors[6],
size = 1,
alpha=.6)+
ggnewscale::new_scale_fill() +
geom_sf(data = recortedatalum,size=0.1,alpha=.6)+
#scale_fill_manual(values = colors2) +
#facet_grid(ANIO,~SEXO) +
#facet_grid(SEXO~ANIO)+
#facet_grid(CLASE~SEXO) +
theme_map() +
coord_sf()
# coord_sf(xlim = xlim, ylim = ylim)
mapamx
#mapamx + theme_void() + # Empty theme without axis lines and texts
# theme(
# panel.background = element_rect(fill = "transparent", colour = NA),
# plot.background = element_rect(fill = "transparent", colour = NA),
# legend.background = element_rect(fill = "transparent", colour = NA),
# legend.box.background = element_rect(fill = "transparent", colour = NA),
# legend.position="bottom",
# legend.direction="horizontal",
# text = element_text(family = "Sans Serif", color = "#c2cae3")
#)
#ggpubr
#mapamx + ggpubr::theme_transparent()
outfile <- str_to_lower(glue("images/{map}/{map}_{pal}Mvialosmmxsps.png"))
width = NA
height = 800
asp <- get_aspect_ratio(tor)
if (is.na(width) && !is.na(height)) {
width <- height * asp
} else if (is.na(height) && !is.na(width)) {
height <- width / asp
}
#"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/portraits/torreon/images/torreon"
#https://stackoverflow.com/questions/57868862/crop-unused-margin-in-a-ggplot-plot
#ggsave(outfile, width = 800, height = 800, units = "mm", dpi = "retina",bg="white",limitsize = FALSE)
ggsave(outfile, width = width, height = height, units = "mm", dpi = "retina",bg="white",limitsize = FALSE)Análisis de puntos y zonas.
#### Densidad de puntos
library(spatstat) # to calculate field of point density
#library(maptools) # to convert to point pattern
#create density map
#puntos<-st_geometry(datacom_sub)
#puntos<-recortedatalum_sf
puntos<-recortedatalum
class(puntos)[1] “sf” “tbl_df” “tbl” “data.frame”
spdfpuntos <- as_Spatial(puntos)
class(spdfpuntos)[1] “SpatialPointsDataFrame” attr(,“package”) [1] “sp”
sf_obj <-st_bbox(puntos)
# Convert the sf object to an owin object
window <- as.owin(sf_obj)
sf_bur_fal_coords <- matrix(unlist(puntos$geometry), ncol = 2, byrow = T)
#Then we use the ppp function to create the ppp (point pattern) object using the information from our matrix and the window that we created.
sSp <- ppp(x = sf_bur_fal_coords[,1], y = sf_bur_fal_coords[,2],
window = window, check = T)
si<-bw.ppl(sSp)
#Dens <- density(sSp, adjust = 0.2)
Dens <- density(sSp, sigma=si)
class(Dens)[1] “im”
plot(Dens)
#contour(density(sSp, adjust = 0.2), nlevels = 4)
contour(density(sSp, sigma=si), nlevels = 4)
#catagorise density
Dsg <- as(Dens, "SpatialGridDataFrame")
#class(Dsg)
#plot(Dsg)
Dim <- as.image.SpatialGridDataFrame(Dsg)
#Dcl <- contourLines(Dim)
#class(Dcl)
library(raster)
r <- raster::raster(Dim)
plot(r)
l <- as.image.SpatialGridDataFrame(as(r, "SpatialGridDataFrame"))
levs <- pretty(range(l$z, na.rm = TRUE),nlevels=12)
cl2SLDF <- function (cL, proj4string = CRS(as.character(NA)))
{
if (length(cL) < 1L)
stop("cL too short")
lev_instances <- unlist(lapply(cL, "[[", "level"))
lev_group <- factor(lev_instances)
df <- data.frame(level = levels(lev_group),
value = unique(lev_instances), ## reasonably confident that contourLines returns these in sort order
stringsAsFactors = TRUE)
clcoords <- lapply(cL, function(x) cbind(x[["x"]], x[["y"]]))
m <- nlevels(lev_group)
res <- vector(mode = "list", length = m)
IDs <- paste("C", 1:m, sep = "_")
row.names(df) <- IDs
for (ilev in seq_len(nlevels(lev_group))) {
idx <- which(lev_group == levels(lev_group)[ilev])
res[[ilev]] <- Lines(lapply(clcoords[idx], Line), IDs[ilev])
}
SL <- SpatialLines(res, proj4string = proj4string)
res <- SpatialLinesDataFrame(SL, data = df)
res
}
cL <- contourLines(l, sort(sample(levs)))
cL <- contourLines(l, levels = pretty(range(l$z, na.rm = TRUE), 10))
plot(sf::st_as_sf(cl2SLDF(cL)))
#text(.1, .2, lab = "1")
#SLDF <- sf::st_as_sf(cl2SLDF(cL))
SLDF <- cl2SLDF(cL)
library(sp)
proj4string(SLDF) <- proj4string(spdfpuntos) # assign correct CRS - for future steps
class(SLDF)[1] “SpatialLinesDataFrame” attr(,“package”) [1] “sp”
plot(SLDF)
SLDF <- SLDF[SLDF@data$level != 0,]
SLDFsf <- sf::st_as_sf(SLDF)
library(sf)
## rgeos version: 0.2-19, (SVN revision 394)
## GEOS runtime version: 3.3.8-CAPI-1.7.8
## Polygon checking: TRUE
Polyclust <- st_polygonize(SLDFsf)
plot(Polyclust)
class(Polyclust)[1] “sf” “data.frame”
gas <- st_area(Polyclust)/10000
polysf<-st_as_sf(tor)
Polyclustsf<-st_as_sf(Polyclust)
mapamx<-ggplot()+
geom_sf(data = polysf, inherit.aes = FALSE,color=colors[6],alpha=.10)+
geom_sf(data = puntos,size=.8,alpha=.3,color="grey")+
#geom_sf(data = Polyclustsf, inherit.aes = FALSE,color=colors[6])+
geom_sf(data=Polyclustsf,aes(fill=level,show.legend=TRUE),color="transparent")+
scale_fill_manual(values=colors2) +
#geom_sf(data = SLDFsf, inherit.aes = FALSE,color=colors[5],size=.2)+
geom_sf(data=viasosmlayer_sub |>
filter(fclass %in% c("primary","primary_link")),
inherit.aes = FALSE,
color = colors[5],
size = 0.3,
alpha=.6)+
geom_sf(data=viasosmlayer_sub |>
filter(fclass %in% c("scondary","tertiary","residential")),
inherit.aes = FALSE,
color = colors[3],
size = 0.3,
alpha=.6)+
geom_sf(data=viasosmlayer_sub |>
filter(fclass %notin% c("primary","primary_link","scondary","tertiary","residential")),
inherit.aes = FALSE,
color = colors[8],
size = 0.2,
alpha=.6)+
geom_sf(data=viasredvial_sub |>
filter(CARRILES %in% c("6","5","4")),
inherit.aes = FALSE,
color = colors[1],
size = .2,
alpha=.6)+
geom_sf(data=viasredvial_sub |>
filter(CARRILES %in% c("3","2","1")),
inherit.aes = FALSE,
color = colors[1],
size = .2,
alpha=.6)+
geom_sf(data=viasredvial_sub |>
filter(CARRILES %notin% c("6","5,","4","3","2","1")),
inherit.aes = FALSE,
color = colors[6],
size =.2,
alpha=.6)+
ggnewscale::new_scale_fill() +
#geom_sf(data = recortedatalum,size=0.1,alpha=.6)+
#scale_fill_manual(values = colors2) +
#facet_grid(ANIO,~SEXO) +
#facet_grid(SEXO~ANIO)+
#facet_grid(CLASE~SEXO) +
theme_map() +
coord_sf()
# coord_sf(xlim = xlim, ylim = ylim)
mapamx
width = NA
height = 800
asp <- get_aspect_ratio(tor)
if (is.na(width) && !is.na(height)) {
width <- height * asp
} else if (is.na(height) && !is.na(width)) {
height <- width / asp
}
#"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/portraits/torreon/images/torreon"
#https://stackoverflow.com/questions/57868862/crop-unused-margin-in-a-ggplot-plot
outfile <- str_to_lower(glue("images/{map}/{map}_{pal}Mdensidadpuntos.png"))
#ggsave(outfile, width = 800, height = 800, units = "mm", dpi = "retina",bg="white",limitsize = FALSE)
ggsave(outfile, width = width, height = height, units = "mm", dpi = "retina",bg="white",limitsize = FALSE)