创建一个用不同颜色区分集群的地图

我正在生成一张带有一些房屋位置的地图。但是，请注意，在 中Points_properties，我有一列表示给定房屋所在的集群。在这种情况下，集群在 1 和 2 之间的范围内。我什至制作了一张区分点颜色的地图，但我的想法是在集群的区域中显示不同的颜色，根据下面的示例图，它代表了我想要的想法。

请注意，在上面的地图中，集群 1 区域为绿色，集群 2 区域为蓝色。

可执行代码如下：

    library(rgdal)
    library(sf)
    library(raster)
    
    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)
points(x = Points_properties$Longitude,
       y= Points_properties$Latitude, col = "black",pch = 16,cex=1)