Create a map that differentiates clusters by different colors-CodePudding

The

First code:

library(rgdal) library(sf) library(raster) library (dplyr) temp <- tempfile() temp2 <- tempfile() download.file("https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2015/UFs/PR/pr_municipios.zip",temp) unzip(zipfile = temp, exdir = temp2) shp <- readOGR(temp2) shp_subset <- shp[shp$NM_MUNICIP == "CASTRO",] #Points Points_properties<-structure(list(Latitude = c(-24.781624, -24.775017, -24.769196, -24.761741, -24.752019, -24.748008, -24.737312, -24.744718, -24.751996, -24.724589, -24.8004, -24.796899, -24.795041, -24.780501, -24.763376, -24.801715, -24.728005, -24.737845, -24.743485, -24.742601, -24.766422, -24.767525, -24.775631, -24.792703, -24.790994, -24.787275, -24.795902, -24.785587, -24.787558, -24.799524), Longitude = c(-49.937369, -49.950576, -49.927608, -49.92762, -49.920608, -49.927707, -49.922095, -49.915438, -49.910843, -49.899478, -49.901775, -49.89364, -49.925657, -49.893193, -49.94081, -49.911967, -49.893358, -49.903904, -49.906435, -49.927951, -49.939603, -49.941541, -49.94455, -49.929797, -49.92141, -49.915141, -49.91042, -49.904772, -49.894034, -49.86651), cluster = c("1", "1", "1", "1", "2", "2", "2", "2", "2", "2", "1", "1", "1", "1", "1", "1", "2", "2", "2", "2", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1")), row.names = c(NA, -30L), class = c("tbl_df", "tbl", "data.frame")) #Generate map marsize <- .2 par(mar = rep(marsize,4)) z <- .4 shp_subset_cropped <- raster::crop(shp_subset, extent(c( (extent(shp_subset)@"xmin"*(1-z) z*min(Points_properties$Longitude)), (extent(shp_subset)@"xmax"*(1-z) z*max(Points_properties$Longitude)), (extent(shp_subset)@"ymin"*(1-z) z*min(Points_properties$Latitude)), (extent(shp_subset)@"ymax"*(1-z) z*max(Points_properties$Latitude)) ))) plot(shp_subset_cropped) pp1 <- filter(Points_properties, cluster==1) pp2 <- filter(Points_properties, cluster==2) points(x = pp1$Longitude, y= pp1$Latitude, col = "red",pch = 16,cex=1) points(x = pp2$Longitude, y= pp2$Latitude, col = "blue",pch = 16,cex=1)

Generated output 1

Second code:

library(rgdal) library(ggplot2) library(sf) temp <- tempfile() temp2 <- tempfile() download.file("https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2015/UFs/PR/pr_municipios.zip",temp) unzip(zipfile = temp, exdir = temp2) shp <- sf::read_sf(temp2) shp_subset <- shp[shp$NM_MUNICIP == "CASTRO",] #Points Points_properties<-structure(list(Latitude = c(-24.781624, -24.775017, -24.769196, -24.761741, -24.752019, -24.748008, -24.737312, -24.744718, -24.751996, -24.724589, -24.8004, -24.796899, -24.795041, -24.780501, -24.763376, -24.801715, -24.728005, -24.737845, -24.743485, -24.742601, -24.766422, -24.767525, -24.775631, -24.792703, -24.790994, -24.787275, -24.795902, -24.785587, -24.787558, -24.799524), Longitude = c(-49.937369, -49.950576, -49.927608, -49.92762, -49.920608, -49.927707, -49.922095, -49.915438, -49.910843, -49.899478, -49.901775, -49.89364, -49.925657, -49.893193, -49.94081, -49.911967, -49.893358, -49.903904, -49.906435, -49.927951, -49.939603, -49.941541, -49.94455, -49.929797, -49.92141, -49.915141, -49.91042, -49.904772, -49.894034, -49.86651), cluster = c("1", "1", "1", "1", "2", "2", "2", "2", "2", "2", "1", "1", "1", "1", "1", "1", "2", "2", "2", "2", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1")), row.names = c(NA, -30L), class = "data.frame") #Generate map ggplot() geom_sf(data = shp_subset, aes()) geom_point(data = Points_properties, aes(x = Longitude, y = Latitude, color = cluster)) coord_sf(xlim = c(min(Points_properties$Longitude)-0.1, max(Points_properties$Longitude) 0.1), ylim = c(min(Points_properties$Latitude)-0.1, max(Points_properties$Latitude) 0.1), expand = FALSE) theme_void()

Generated output 2

Map generated by Phil code

CodePudding user response：

This is imperfect because the fill extends beyond the boundaries of the shapefile, but it might be useful nonetheless. The trick is to use a voronoi algorithm using ggraph::geom_node_voronoi().

library(ggraph) ggplot() geom_sf(data = shp_subset) geom_point(data = Points_properties, aes(x = Longitude, y = Latitude, color = cluster)) geom_node_voronoi(data = Points_properties, bound = c(-50.24633, -49.49878 , -25.07238, -24.541), aes(x = Longitude, y = Latitude, fill = cluster), alpha = 0.2) coord_sf(xlim = c(min(Points_properties$Longitude)-0.1, max(Points_properties$Longitude) 0.1), ylim = c(min(Points_properties$Latitude)-0.1, max(Points_properties$Latitude) 0.1), expand = FALSE) theme_void()

CodePudding user response：

Using terra (the replacement of raster)

Your example data

library(terra) mun <- vect(file.path(temp2, "41MUE250GC_SIR.shp")) mun <- mun[mun$NM_MUNICIP == "CASTRO",] pts <- vect(data.frame(Points_properties), c("Longitude", "Latitude"))

Create a Voronoi diagram for the extent of the municipality, aggregate it, and remove areas outside the municipality (I am guessing that is what you want)

z <- voronoi(pts, mun) z <- aggregate(z, "cluster") z <- crop(z, mun)

Now plot

plot(z, "cluster", col=c("blue", "green")) points(pts, col="red", cex=.75)