如何在R语言spplot绘制的地图中为指定区域添加阴影多边形？

我有一份统计巴黎各街区房地产饱和度百分比的地图，指标规则为房地产实际观测价格/法律规定的最高限价。我需要为数据集中观测数不足5的街区添加阴影填充效果。我检索了相关资料但未找到可行实现方案，欢迎各位提供正确的实现思路，非常感谢。

现有实现代码 #

library(sp)
library(sf)
library(rgdal)
library(RColorBrewer)
library(raster)
library(classInt)
library(cartography)

#Importation
setwd("path")
shp <- readOGR(dsn="path/to/file",layer="l_qu_paris")
#Breaks
q10 <- classIntervals(map$saturation2, n=7, style="fixed",
               fixedBreaks=c(45,69.999999, 79.9999999, 89.9999999, 99.9999999
                             ,109.99999999, 120))
#Colors
my.palette <- colors()[c(73,26,128,10,652,92)]
#Map
##Scale
scale.parameter = 1.1
xshift = 0
yshift = 0
original.bbox = shp@bbox
edges = original.bbox
edges[1, ] <- (edges[1, ] - mean(edges[1, ])) * scale.parameter + mean(edges[1, ]) + xshift
edges[2, ] <- (edges[2, ] - mean(edges[2, ])) * scale.parameter + mean(edges[2, ]) + yshift

#Saturation
idx <- match(shp$l_qu, map$l_qu)
is.na(idx)
concordance <- map[idx, "saturation2"]
shp$saturation2 <- concordance
spplot(shp, "saturation2",col.regions=my.palette, 
       col = "black", lwd= 1, at = q10$brks,
       main=list(label="% de saturation des meublés 1 pièce",cex=1.2,fontfamily="serif"),
       xlim = edges[1, ], ylim = edges[2, ])
grid.text("Saturation moyenne (en%)", x=unit(0.95, "npc"), y=unit(0.50, "npc"), rot=90)

所用数据集 #

map <- structure(list(l_qu = c("Amérique", "Archives", "Arsenal", "Arts-et-Métiers", 
"Auteuil", "Batignolles", "Bel-Air", "Belleville", "Bercy", "Bonne-Nouvelle", 
"Chaillot", "Champs-Elysées", "Charonne", "Chaussée-d'Antin", 
"Clignancourt", "Combat", "Croulebarbe", "Ecole-Militaire", "Enfants-Rouges", 
"Epinettes", "Europe", "Faubourg-du-Roule", "Faubourg-Montmartre", 
"Folie-Méricourt", "Gaillon", "Gare", "Goutte-d'Or", "Grandes-Carrières", 
"Grenelle", "Gros-Caillou", "Halles", "Hôpital-Saint-Louis", 
"Invalides", "Jardin-des-Plantes", "Javel", "La Chapelle", "Madeleine", 
"Mail", "Maison-Blanche", "Monnaie", "Montparnasse", "Muette", 
"Necker", "Notre-Dame", "Notre-Dame-des-Champs", "Odéon", "Palais-Royal", 
"Parc-de-Montsouris", "Père-Lachaise", "Petit-Montrouge", "Picpus", 
"Place-Vendôme", "Plaine de Monceaux", "Plaisance", "Pont-de-Flandre", 
"Porte-Dauphine", "Porte-Saint-Denis", "Porte-Saint-Martin", 
"Quinze-Vingts", "Rochechouart", "Roquette", "Saint-Ambroise", 
"Saint-Fargeau", "Saint-Germain-des-Prés", "Saint-Gervais", 
"Saint-Lambert", "Saint-Merri", "Saint-Thomas-d'Aquin", "Saint-Victor", 
"Saint-Vincent-de-Paul", "Sainte-Avoie", "Sainte-Marguerite", 
"Saint-Georges", "Salpêtrière", "Sorbonne", "Saint-Germain-l'Auxerrois", 
"Ternes", "Val-de-Grâce", "Villette", "Vivienne", "Total"), 
    saturation2 = c(98.188951329533, 85.4938271604938, 83.8463463463464, 
    90.1460755525873, 98.1726527090667, 90.2186740262059, 92.8743271072797, 
    72.8549079897508, 99.2356140350877, 90.1234567901235, 114.057904044022, 
    NA, 87.2208980972528, 91.2562612612613, 97.9518951016991, 
    86.2770900920801, 91.0239726151895, 92.8305400372439, 88.6514719848053, 
    73.876877752942, 108.693318725755, 67.3263578578579, 85.8735259484408, 
    89.2100224414912, 92, 90.6120989320281, 85.8446948520848, 
    91.4165103088783, 97.2760978594495, 93.60892313074, 102.471730530348, 
    95.9062868379746, 96, 92.5484278273071, 95.0066946433545, 
    85.8187074829932, 101.139150713213, 92.1272297297297, 93.0625144594594, 
    61.8074324324324, 100.173302938197, 99.720856146949, 84.8732544128823, 
    84.1911355800245, 85.1122672253259, 91.8422003734504, NA, 
    94.612349767814, 83.2363741480137, 87.0403187718064, 92.0886931496388, 
    77, 110.943302180685, 100.73486307088, 66.3899425287356, 
    96.2527514568292, 95.7430893746874, 87.9028997984617, 48, 
    85.5630809345015, 92.7010730078939, 82.075822827797, 83.1727736726875, 
    76.2162162162162, 104.534662867996, 98.3510353194912, 78.3333333333333, 
    103.169134078212, 80.8779605984059, 92.63515704154, 62, 90.3902768982325, 
    94.1391771653151, 94.8669917042241, 94.4825319797959, 95.4279279279279, 
    98.2238673533848, 94.0602977590835, 87.5105365473892, 102, 
    92.5123935729199), numobs = c(6, 4, 4, 6, 36, 15, 4, 4, 3, 
    2, 16, NA, 36, 3, 32, 9, 22, 13, 11, 6, 31, 5, 15, 14, 4, 
    22, 3, 64, 29, 58, 7, 18, 4, 13, 23, 2, 8, 4, 47, 12, 16, 
    49, 50, 9, 33, 26, NA, 15, 10, 10, 23, 2, 13, 15, 2, 12, 
    8, 31, 1, 17, 22, 42, 7, 3, 4, 74, 4, 7, 13, 6, 2, 23, 18, 
    16, 17, 1, 24, 44, 8, 4, 1290)), row.names = c(NA, -81L), class = c("tbl_df", 
"tbl", "data.frame"))

内容的提问来源于stack exchange，提问作者MonsieurDjuna